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ABSTRACT 


'''This  report  contains  a  description  of  the  work  done  in  the  first  phase  of  a 
contract  awarded  by  the  Defence  Research  Establishment  Atlantic  (DREA)  for  the 
provision  of  an  in-house  random  response  analysis  capability.  A  concise 
description  of  the  theoretical  basis  of  the  finite  element  random  vibration  analysis 
methodology  is  given.  A  description  of  available  procedures  for  computing  the 
assembled  cross-spectral  density  matrix  for  a  finite  element  model  of  structures 
excited  by  distributed  random  loads  is  also  presented.  Special  emphasis  is  given 
to  marine  structures  subjected  to  random  pressure  loads  induced  by  random 
ocean  waves.  Computer  programs  for  performing  random  vibration  analysis 
developed  during  the  course  of  this  work  have  been  incorporated  into  the  VAST 
finite  element  program  as  a  new  module  called  RANVIB.  A  suite  of  computer 
programs  for  the  generation  of  random  loads  for  ship  structures  were  also 
developed.  Descriptions  of  all  programs  together  with  operational  procedures  and 
input  data  preparation  are  provided.  Example  problems  illustrating  the  operation 
of  RANVIB  are  given.  The  report  also  gives  some  remarks  about  the  capabilities 
developed  and  recommendations  for  future  work. 


RESUME 


Le  present  rapport  decrit  les  travaux  effectues  au  cours  de  la  premiere 
phase  dun  contrat  attribue  par  le  Centre  de  recherche  pour  la  defense  Atlantique 
(CRDA)  pour  la  creation  d’un  service  interne  d’analyse  des  reponses  aleatoires. 
Une  description  concise  de  la  base  theorique  de  la  methode  d’analyse  des 
vibrations  aleatoires  par  elements  finis  est  donn<§e.  Une  description  des  methodes 
disponibles  pour  calculer  la  matrice  des  densites  trans  spectrales  elaboree  pour 
un  modele  aux  elements  finis  de  structures  excitees  par  des  charges  aleatoires  est 
aussi  presentee.  L'accent  est  mis  sur  les  ouvrages  marins  soumis  a  des 
pressions  aleatoires  exercees  par  des  vagues  oc4aniques  aleatoires.  Les  logiciels 
d'analyse  des  vibrations  aleatoires,  elabor^s  au  cours  de  ces  travaux,  ont  ete 
incorpores  dans  le  programme  aux  elements  finis  VAST  sous  fourme  d’un 
nouveau  module  appeie  RANVIB.  Une  suite  de  logiciels  generant  des  charges 
aleatoires  pour  des  constructions  navales  ont  aussi  ete  eiabores.  Des  descriptions 
de  tous  les  programmes  ainsi  que  des  modes  operatoires  et  du  traitement  des 
donnees  d'entree  sont  fournies.  Des  exemples  illustrant  le  fonctionnement  de 
RANVIB  sont  donnes.  Le  rapport  renferme  aussi  quelques  remarques  sur  les 
competences  acquises  et  des  recommandations  de  travaux  futurs. 
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l.o  INTRODUCTION 


The  fact  that  uncertainty  is  inherent  in  virtually  all  natural  and 
man-made  phenomena  is  widely  accepted.  The  use  of  deterministic  concepts 
such  as  the  factor  of  safety  to  account  for  uncertainties  associated  with 
some  problems  is  inadequate  for  the  cases  in  which  the  level  of 
uncertainty  is  high  or  where  a  high  degree  of  precision  is  needed  to 
describe  system  behaviour.  The  probabilistic  approach  to  system 
uncertainties  provides  a  very  good  framework  for  dealing  with  this  class 

of  problems  because  it  is  based  on  a  rational  and  realistic  methodology 
for  dealing  with  such  uncertainties. 


In  the  field  of  structural  dynamics,  in  particular,  considerable 
efforts  have  been  directed  toward  the  computation  of  the  response  of 
engineering  structures  subjected  to  random  excitations  over  the  years. 
This  is  because  many  modern  structures  and  their  components  must  be 
designed  to  withstand  a  variety  of  excitations  that  are  best  described  as 
nondeterministic.  Such  structures  include  aircrafts  flying  in  gusty 
atmospheres;  high-rise  buildings,  suspension  bridges,  and  suspended  roofs 
subjected  to  wind  and  earthquake  loads;  aerospace  structures  and  missiles 
subjected  to  jet  noise  excitations;  vehicles  travelling  on  rough  tracks; 
ships,  submarines  and  offshore  structures  subjected  to  the  hostile 
environment  of  the  confused  sea;  aircrafts  landing  on  imperfect  runways; 
and  heat  exchanger  tubes  subject  to  turbulent  flows. 


Primarily,  it  has  been  the  applications  in  aerospace  engineering 

which  have  sped  the  development  of  the  probabilistic  techniques  in  the 

area  of  structural  dynamics.  However,  these  techniques  are  being 

increasingly  accepted  and  applied  for  the  design  and  analysis  of  the 

integrity  of  other  engineering  structures,  including  ship  and  offshore 
structures. 
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Ocean  waves,  which  constitute  major  loadings  for  ships  and  offshore 
platforms,  are  inherently  random  in  terms  of  both  temporal  and  spatial 
characteristics.  As  such,  stochastic  modelling  plays  an  important  role 
in  the  design,  operation  and  service  of  these  marine  structures.  The 
probabilistic  characterizations  and  determinations  of  their  responses 
have  attracted  the  attention  of  scientists  and  engineers.  A  number  of 
interesting  studies  addressing  this  subject  area  have  been  published,  see 
for  example  [1,2,3]. 

In  realization  of  the  need  to  adopt  a  probabilistic  approach  for  the 
design  and  structural  integrity  assessment  of  Canadian  naval  ships,  the 
Defence  Research  Establishment  Atlantic  (DREA)  has  decided  to  provide  an 
in-house  capability  to  perform  the  required  computations.  Three 
important  components  are  required  for  this  capability.  First,  an 
adequate  description  of  the  random  wave  loads  imparted  on  the  ship 
structures  by  the  random  sea  surface  elevation  is  needed.  Secondly,  the 
capability  to  determine  the  random  vibration  response  to  the  random 
excitations  in  a  finite  element  setting  must  also  be  available.  Finally, 
provision  must  be  made  for  the  utilization  of  the  response  statistics 
available  in  the  second  capability  for  the  estimation  of  the  reliability 
of  the  ship  structure.  This  third  and  final  component  is  the  ultimate 
goal  of  any  probabilistic  approach  to  design  or  the  assessment  of  the 
integrity  of  an  engineering  structure. 

A  two-phase  contract  was  recently  awarded  to  Martec  Limited  to  meet 
some  of  the  objectives  described  in  the  preceding  paragraph.  In  the 
first  phase,  a  methodology  for  describing  the  random  pressure  loads 
imparted  on  ship  structures  travelling  in  rough  seas  with  known  wave 
spectra  is  to  be  provided.  The  random  loading  is  to  be  described  in  a 
form  suitable  for  input  into  the  finite  element  program  that  will  be  used 
for  the  computation  of  the  random  response  quantities.  Secondly,  a 
random  vibration  analysis  capability  is  to  be  provided  in  the  finite 
element  program  VAST.  In  the  second  phase  of  the  work,  a  state-of-the- 
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art  literature  survey  is  to  be  carried  out  for  the  purpose  of  identifying 
suitable  probabilistic  methodologies  to  be  adopted  for  the  accurate 
(fatigue)  reliability  assessment  of  naval  ships. 

The  present  document  represents  the  final  report  of  the  first  phase 
of  the  contract  work.  In  Chapter  2,  a  brief  introduction  to  the  theory 
of  Random  Vibration  is  given  especially  the  part  that  Is  of  relevance  to 
the  one  implemented  in  the  current  work.  In  Chapter  3,  a  description  of 
the  random  load  generation  is  presented.  A  detailed  description  of  the 
implementation  of  the  random  vibration  analysis  module  in  VAST  called 
RANVIB  is  presented  in  Chapter  4.  A  description  of  the  set  of  computer 
programs  developed  for  the  computation  of  the  random  wave  loads  imposed 
on  ship  hulls  is  also  given.  The  provision  of  a  graphics  support  for  the 
RANVIB  module  is  given  in  Chapter  5.  Example  problems  used  for  the 
verification  of  program  RANVIB  and  illustrating  its  use  are  presented  in 
Chapter  6.  Chapter  7  closes  the  report  by  giving  a  few  remarks  about  the 
work  done  during  the  course  of  this  work  and  some  recoimiendations  for 
additional  work  that  would  enhance  the  versatility  of  the  capabilities 
that  have  been  provided. 
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2.0  RANDOM  VIBRATION  THEORY 


The  subject  of  random  vibrations  is  concerned  with  the  application 
of  probabilistic  methods  for  the  determination  of  the  response  of  certain 
or  uncertain  systems  subjected  to  dynamic  loads  that  are  best  described 
as  nondeterministic.  As  a  technical  discipline,  therefore,  random  vibra¬ 
tion  is  a  combination  of  structural  dynamics  and  probability  theory. 

A  dynamic  excitation  representable  as  a  function  of  time  is  said  to 
be  nondeterministic  or  random  if  the  function  cannot  be  given  any  explic¬ 
it  certain  description.  In  other  words  the  variation  of  the  value  of  the 
function  with  time  is  so  irregular  and  uncertain  that  it  can  only  be 
described  in  a  statistical  or  averaged  sense  and  so  the  excitation  is 
said  to  constitute  a  random  process. 


If,  in  addition  to  randomness  in  time,  there  is  also  spatial  random¬ 
ness,  the  excitation  is  referred  to  as  a  random  field.  Thus,  the  key 
concepts  from  probability  theory  are  those  of  random  processes  and  random 

fields  which  are  used  as  stochastic  models  for  excitation  and  response 
histories. 

2.1  Classification  of  Random  Processes 

Based  on  a  statistical  regularity  classification  scheme  [4],  time- 
parametered  random  processes  are  generally  grouped  into  two  main  categor¬ 
ies:  stationary  and  nonstationary.  Analogously,  space-parametered  random 
processes  are  classified  as  homogeneous  or  nonhomogeneous. 

A  stationary  process  is  one  in  which  the  statistical  behaviour  does 
not  vary  with  time.  Its  complete  probability  structure,  which  for  a 
stochastic  process  is  defined  by  a  system  of  probability  distribution 
functions,  is  independent  of  a  shift  in  the  time  origin  across  the 
ensemble.  Thus,  the  properties  of  this  process  can  hypothetically  be 
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described  at  any  instant  of  time  by  computing  average  values  over  the 
collection  of  sample  functions  which  describe  the  random  process  and  they 
will  be  found  to  be  time-invariant.  For  an  excitation  that  can  be  so 
classified,  its  duration  is  usually  relatively  long  compared  with  the 
period  of  the  lowest  frequency  spectral  components  of  the  excited  system 
[5].  Similarly,  a  random  field  is  said  to  be  homogeneous  with  respect  to 
a  particular  spatial  coordinate  if  its  probability  distributions  are 
invariant  with  respect  to  translations  of  the  origin  along  the  axis  of 
that  coordinate. 

On  the  other  hand,  a  nonstationary  (nonhomogeneous)  process  is  one 
whose  probability  distribution  functions  depend  upon  the  values  of  time 
(space)  parameters  explicitly.  In  other  words,  its  statistical  behaviour 
depends  upon  the  absolute  origin  of  time  (space).  Furthermore,  a  non¬ 
stationary  (nonhomogeneous)  excitation  is  such  that  its  duration  is  not 
much  longer  than  the  fundamental  time  (spatial)  period  of  the  system. 
Nonstationary  random  processes  are  in  general  much  more  difficult  to 
handle  than  their  stationary  counterparts:  measurement  or  description, 
response  prediction,  and  drawing  conclusions  about  reliability.  The 
majority  of  disturbances  encountered  in  real-life  applications  are 
actually  nonstationary  by  definition.  This  is  because  all  engineering 
random  process  must  have  a  beginning  and  ending.  In  many  cases,  however, 
it  is  acceptable  to  assume  a  kind  of  uniformity  or  regularity  for  the 
random  field  or  random  process  describing  such  disturbances  for  the 
lifetime  or  piecewise  portions  thereof.  Fortunately,  the  case  of 
interest  in  this  work  falls  into  the  category  of  problems  that  can 
reasonably  be  assumed  to  be  stationary  in  time  and  homogeneous  in  space. 
As  such,  all  further  discussions  of  the  theory  of  random  vibration  will 
focus  on  stationary  excitations  and  response  computations. 

A  stationary  random  process  may  be  further  classified  as  ergodic  or 
non-ergodic.  An  ergodic  process  is  one  in  which,  in  addition  to  all  the 
ensemble  averages  being  stationary  with  respect  to  a  change  of  the  time 
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scale,  the  averages  taken  along  any  single  sample  are  the  same  as  the 
ensemble  averages.  In  practical  terms,  each  sample  function  of  an 
ergodic  process  is  completely  representative  of  the  ensemble  that 
constitutes  the  random  process.  For  a  non-ergodic  stationary  process, 
the  averages  over  a  single  realization  of  the  random  phenomenon  cannot  be 
used  to  describe  the  ensemble  average. 

2.2  Introduction  to  Random  Vibration  Terminology 

For  completeness  of  this  report,  the  definitions  of  a  few  Random 
Vibration  terminologies  are  given  in  what  follows. 

Let  X(t)  represent  a  random  process  that  is  time  stationary.  The 
probability  distribution  function  of  X(t),  denoted  Fx(x),  is  the  proba¬ 
bility  that  X  assumes  a  value  less  than  x,  i.e. 

Fx(x)  =  P(X  <  x)  (2.1) 

In  practice  it  is  more  convenient  to  deal  with  the  derivatives  of  the 
probability  distribution  functions,  called  the  probability  density 
functions,  px(x),  defined  as: 


Px(x) 


dFy(x) 

dx 


(2.2a) 


Thus,  for  a  process  X(t),  p(x)dx  is  the  fraction  of  the  total  elapsed 
time  for  which  X(t)  lies  in  the  band  x  to  x+dx.  Alternatively,  equation 
(2.2a)  may  be  expressed  as 

x 

P(X<x)  =  /  p(e)de 


(2.2b) 
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The  expected  value  of  a  random  variable  x  is  defined  as 


co 

E[X]  =  /  xpx(x}dx 


(2.3) 


An  expected  value  is  synonymously  called  an  ensemble  average  (or  mean), 
or  statistical  average  (or  mean),  or  mathematical  expectation. 

The  mean  square  of  X  is  given  by 


E[X2]  =  /  x2px(x)dx 


(2.4) 


and  the  variance,  o2,  is  given  by 

o2  =  E[(X-E[XH)2]  (2.5a) 
or  o2  =  E[X2]  -  (E[X])2  (2.5b) 

The  mean  and  the  mean  square  of  a  random  variable  or  a  random 
process  depend  on  the  first-order  probability  distribution  as  can  be  seen 
in  equations  (2.3)  and  (2.4).  These  are  the  most  basic  and  common 
statistical  averages  used  in  the  description  of  random  processes.  For  a 
zero-mean  process,  the  variance  is  the  same  as  the  mean  square  value 
which  gives  a  rudimentary  description  of  the  intensity  of  the  random 
process. 

While  the  probability  density  function  p(x)  gives  an  amplitude 
domain  description  of  the  process,  the  autocorrelation  and  power  spectral 
density  functions  furnish  similar  information  in  the  time  and  frequency 
domains,  respectively.  If  f(xl5tj)  is  a  random  field,  the  space-time  cor¬ 
relation  function  is  the  ensemble  average  of  the  product  f(xx ,t,)f (x2 ,t2) . 
In  general,  this  is  a  function  of  the  two  locations  and  the  two  times.  If 
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the  field  is  stationary  then  the  correlation  is  merely  dependent  on  the 
difference  t  =  t2-ti  of  the  two  times.  Similarly,  for  a  homogeneous 
process,  the  correlation  depends  only  on  the  spatial  distance  x2-x!. 

For  a  time-stationary  process,  the  autocorrelation  function,  R  (t) 
is  defined  as  XX 

Rxx(t)  =  E[x(t)x(t+t)]  (2.6a) 

while  the  cross-correlation  function  of  two  processes  x(t)  and  y(t)  is 
defined  as 

Rxy(t)  =  E[x(t)y(t+x)]  (2.6b) 

The  Fourier  transform  of  the  correlation  function  of  a  stationary  process 
is  called  the  cross-spectral  density  function  and  is  mathematically 
expressed  as 


Sxy(u)  =  h  RXy(x)exP(-i(,)T)dT  (2.7) 

Thus,  when  the  cross  spectral  density  function  is  known,  the  cross 
correlation  function  may  be  readily  obtained  using  the  expression: 


Rxy(t)  =  J  SXy (to)exp(i tux ) dto  (2.8) 

Of  particular  practical  significance  is  that  the  value  of  the 
autocorrelation  function  evaluated  at  x=0  gives  the  mean  square  value, 
i.e.: 


E[X2]  =  Rxx(x=0)  =  /  Sxx(co)dw 


(2.9) 
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Equation  (2.9)  represents  the  area  under  the  graph  of  power  spectrum 
versus  frequency. 

The  two  Fourier  transform  pairs  in  (2.7)  and  (2.8)  are  the  well- 
known  Wiener-Khintchine  relations  [6,7].  The  importance  of  the  correla¬ 
tion  and  spectral  density  functions  as  statistical  descriptors  of  random 
processes  stem  from  the  fact  that  they  provide  average  amplitude  and 
frequency  information  about  the  sample  histories.  Furthermore,  they  can 
be  measured  with  available  data-processing  techniques;  they  are  mathe¬ 
matically  closed  with  respect  to  linear  time-invariant  operations  in  the 
sense  that  if  these  statistics  are  known  for  the  excitation,  then  it  is 
possible  to  obtain  the  corresponding  statistics  for  the  response  of  a 
linear  time-invariant  dynamic  system;  and,  finally,  they  often  provide 
adequate  information  about  the  response  for  making  vital  engineering 
decisions  concerning  the  severity  of  the  vibration,  extent  of  deforma¬ 
tion,  and  the  reliability  or  safety  of  the  system  [8].  Note  that 
definitions  of  (2.7)  and  (2.8)  vary  in  that  the  factor  l/2u  or  other 
factors  may  be  used  as  a  multiplier  in  Equation  (2.8). 

The  correlation  function  is  a  measure  of  the  dependence  between  two 
random  processes.  It  also  gives  information  about  the  frequencies 
present  in  a  random  process  indirectly  as  can  be  seen  by  examination  of 
the  Wiener-Khintchine  relations  given  above.  The  time  history  x(t)  of  a 
sample  function  of  a  naturally  occurring  random  process  is  not  periodic 
and  so  it  cannot  be  represented  by  a  discrete  Fourier  series.  Further¬ 
more,  for  a  stationary  process,  x(t)  goes  on  forever  and  the  integra- 
bility  condition 

co 

/  I x(t) I dt  <  -  (2.10) 

_oo 

is  not  satisfied.  This  is  why  the  frequency  content  of  a  signal  is 
studied  either  from  its  autocorrelation  function  which  is  amenable  to 
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classical  Fourier  methods  or  a  truncated  version  xT(t)  of  the  random 
process.  Until  recently,  the  first  approach  just  mentioned  was  the  one 
usually  adopted  for  estimating  the  spectrum  of  a  discrete  time  series. 
However,  the  second  approach  is  now  heavily  favoured  in  modern  signal 
processing  and  analysis  since  the  development  of  the  Fast  Fourier  Trans¬ 
form  (FFT)  which  has  revolutionized  the  methodology  for  performing  dis¬ 
crete  Fourier  transforms.  In  this  latter  approach,  rather  than  estimate 
the  spectrum  by  first  determining  the  autocorrelation  function  and  then 
calculating  the  Fourier  transform,  the  spectra  are  more  efficiently  and 
more  accurately  computed  directly  from  the  original  time  signal.  We 
shall  see  later  that  a  knowledge  of  the  spectrum  of  a  random  excitation 
is  very  vital  in  the  application  of  frequency  domain  methods  for  the 
calculation  of  random  responses  to  such  excitations. 

2,3  Stationary  Random  Vibration  Analysis  of  Linear  Systems 

Random  vibration  analysis  involves  the  use  of  probabilistic  methods 
to  perform  structural  dynamic  analysis.  The  probabilistic  description  of 
the  random  excitation  is  used  to  determine  the  statistical  properties  of 
the  random  response.  In  general,  the  response  of  a  linear  time-invariant 
system  to  a  stationary  excitation  is  also  stationary.  The  only  exception 
to  this  rule  is  if  the  system  has  no  damping  or  is  unstable  [8].  In  the 
application  which  is  of  interest  in  this  work,  it  is  very  safe  to  consid¬ 
er  our  random  response  to  be  stationary  because  all  real-life  structures 

possess  some  degree  of  damping  and  we,  of  course,  assume  our  structure  to 
be  stable. 

In  conventional  vibration  or  structural  dynamics  analysis,  there  are 
two  main  techniques  for  the  solution  of  dynamic  problems.  In  the  first 
technique,  the  equations  of  motion  are  directly  integrated  using  either 
implicit  or  explicit  time-stepping  numerical  integration  schemes  to 
determine  the  dynamic  response.  This  technique  is  popularly  referred  to 
as  the  time  domain  method.  In  the  second  approach,  the  power  of  the 
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Fourier  series  for  period  functions  or  Fourier  transforms  for  nonperiodic 
functions  is  applied.  The  beauty  of  this  approach  for  linear  systems 
lies  in  the  applicability  of  the  principle  of  superposition  whereby  the 
response  of  the  dynamic  system  is  obtained  as  the  summation  of  the 
responses  of  each  of  the  terms  of  the  Fourier  series  expansion  or  the 
discrete  Fourier  transform  of  the  excitation  function.  This  latter  tech¬ 
nique  known  as  the  frequency  domain  approach  is  particularly  attractive 
for  linear  systems  because  the  numerical  integration  involved  in  the 
direct  method  is  replaced  by  finite  sums  by  using  the  powerful  FFT.  The 
frequency  domain  approach  may  be  further  subdivided  into  the  direct 
frequency  response  method  and  the  modal  frequency  response  method. 
Linear  dynamic  analysis  using  modal  superposition  (or  the  normal -mode 
method)  is  computationally  inexpensive  apart  from  providing  useful 
insight  into  the  dynamic  behaviour  of  the  system.  Details  of  the  above 
methods  can  be  found  in  any  classical  textbook  on  vibrations  or  structur¬ 
al  dynamics. 

Similar  to  conventional  vibration  analysis,  the  computation  of  the 
stationary  random  response  of  a  system  subjected  to  stationary  random 
excitation  can  proceed  via  a  time-domain  approach  or  a  frequency-domain 
approach.  The  important  difference  is  that  the  equations  of  motion  (or 
their  Fourier  transforms  in  the  case  of  the  frequency-domain  approach) 
are  dealt  with  in  a  probabilistic  sense  through  the  statistical  proper¬ 
ties  of  the  excitation.  In  the  time-domain  case,  the  auto-  and  cross¬ 
correlation  functions  represent  the  description  of  the  applied  forces  in 
time  while  the  auto-  and  cross-spectral  densities  are  the  frequency- 
domain  descriptors  of  the  statistical  properties  of  the  forces.  As  was 
mentioned  earlier  on,  the  modern  use  of  the  powerful  FFT  technique  has 
made  it  possible  to  directly  obtain  the  spectral  densities  of  random 
processes.  As  such,  the  application  of  a  frequency-domain  approach  to 
the  random  response  analysis  of  linear  structures  subjected  to  stationary 
loads  is  very  attractive. 
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2*3.1  Response  of  a  Single  Degree-of-Freedom  System 

Consider  a  single  degree-of-freedom  (SDF)  system  comprising  of  mass 
m,  viscous  damping  coefficient  c,  stiffness  k,  and  whose  motion  is 
governed  by  the  equation: 


mx  +  cx  +kx  =  F(t)  (2.11) 

where  F(t)  is  the  stationary  random  force  applied  to  the  mass,  x=x(t)  is 
the  corresponding  stationary  random  displacement  response,  x(t)  is  the 

velocity  and  x(t)  is  the  acceleration.  Let  F(ioj)  denote  the  Fourier 
transform  of  F(t)  and  X(io>)  denote  the  Fourier  transform  of  x(t).  Taking 
the  Fourier  transform  of  both  sides  of  Equation  (2.11)  gives: 


(-u2m  +  iwc  +  k)X(iio)  =  F(ico) 

or  X ( io>)  =  (k  -  to2 m  +  iwc)-1  F(iw) 

or  X(iw)  =  H(iw)F(iw) 
where 


(2.12a) 

(2.12b) 

(2.12c) 


H(iw)  =  (k  -  w2m  +  iwc)-1  (2.12d) 

is  the  well-known  complex  frequency  response  function  or  receptance  of 
the  system. 


If  the  random  function  F(t)  is  truncated  so  that  it  becomes  zero 
outside  the  interval  (-T,+T)  to  give  a  new  function  Fj(t),  then  the 
power  spectral  density  of  F(t)  can  be  expressed  in  the  form  [9] 


SFF(u)  -  |i!f  l[FTdo))]l2  (2.13a) 

sFf(w)  =  T  [PjCiw) ]  •  [F*  (iw)]  (2.13b) 

where  F*(ioj)  is  the  complex  conjugate  of  FT(iw).  Henceforth,  FT(iu>)  and 
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XT(icu)  will  be  used  interchangeably  with  F(iu>)  and  X( ico)  respectively. 


In  other  words.  Equation  (2.13b)  may  also  be  written  as: 

SFF(<0)  -  f  [F(1»)][F*(1w)]  (2.14) 

In  a  similar  manner,  the  power  spectrum  of  x(t)  may  be  expressed  as 

Sxx(«)  =  f  [X(ico) ][X*(ico) ]  (2.15) 

Substituting  the  expression  for  X(iw)  (equation  (2.12c))  into  Equation 
(2.17)  gives 

SXX(W)  =  11-  t  [H(1«)F(1«)][H(1w)F(1w)]*  (2.16a) 

or  Sxx(w)  =  H(iw)  •  {f  [F(io>)  ][ F* ( i cu )  ]}  •  H*(iw)  (2.16b) 

Invoking  the  definition  of  SFF(co)  as  given  by  equation  (2.14)  in  equa¬ 
tion  (2.16b)  gives 


Sxx(u)  =  H(ia>)  •  SpF(w)  •  H*(i<o)  (2.17a) 
or  sxx(w)  =  | H( ito)  | 2  •  SFF(u)  (2.17b) 

Equation  (2.17a)  or  (2.17b)  gives  the  all-important  fundamental 
relationship  between  the  power  spectral  density  of  the  forcing  function 
and  that  of  the  corresponding  response  for  a  single  input-single  output 
system.  The  implication  of  this  relationship  is  that  the  statistical 
properties  of  the  response  can  be  determined  from  those  of  the  input 
provided  the  system  frequency  response  function  is  known.  The  mean  value 
can  also  be  shown  to  be  given  by  [8]: 


E[X(t)3  =  H(0)E[F(t)] 


(2.18) 
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It  is  pertinent  to  point  out  that  the  relations  (2.17)  and  (2.18)  are 
independent  of  the  probability  distributions  involved. 

2.3.2  Response  of  Multi  Deqree-of-Freedom  Systems 

The  results  derived  for  the  single  degree-of-freedom  system  in  the 
preceding  subsection  extend  formally  to  multiple  input-multiple  output 
time-invariant  linear  systems  and  so  can  be  readily  applied  to  discrete 
multi  degree-of-freedom  systems  or  continuous  structures  discretized  by  a 
numerical  technique  such  as  the  finite  element  method. 

The  governing  equation  of  motion  of  a  continuous  structure  discret¬ 
ized  by  the  finite  element  method  is  of  the  form: 


MX  +  CX  +  KX  =  F(t)  (2.21) 

In  Equation  (2.21),  M  is  the  assembled  mass  matrix,  C  is  the  damping 
matrix,  K  is  the  assembled  stiffness  matrix,  X(t)  is  the  vector  displace¬ 
ment  response  to  the  vector  excitation  F(t),~X(t)  is  the  random  velocity 
vector,  and  X(t)  is  the  random  acceleration  vector.  For  our  present 
interest,  it  is  assumed  that  the  system  parameters  M,  C,  and  K  are  known 
deterministic  quantities  thereby  leaving  F(t)  the  only  random  input  for 
which  the  random  output  X ( t )  is  sought.  For  a  structure  discretized  in 
this  mannner  with  a  total  of  NS  global  degrees  of  freedom,  M,  C,  and  K 
are  (NSxNS)  square  matrices  while  both  X  and  F  are  vectors  with  NS  com¬ 
ponents. 

It  can  be  shown,  in  a  manner  similar  to  the  single  degree-of-freedom 
case,  that  the  cross-spectral  density  of  the  vector  random  displacement 

reponse  is  related  to  that  of  the  vector  random  force  excitation  via  the 
expression 


(2.22) 


[Sxx(u>)]  =  [H(iw)][SpF  (w)][H(ifo)]T* 

where,  now,  [Spp((o)]  and  [Sxx(<»)]  are  matrices  of  excitation  and 
response  spectra  including  auto-  and  cross-spectral  density  functions, 
and  [H(icu)]  is  the  complex  frequency  response  function  matrix 
corresponding  to  Equation  (2.21).  On  taking  the  Fourier  transform  of 
both  sides  of  Equation  (2.21),  it  is  easy  to  see  that  [H(iu>)]  is  given  by 

CH(iw)]  =  [K  -  (02M  +  iuiC]-1  (2.23) 

The  procedure  for  computing  [Sxx(w)]  by  using  Equations  (2.22)  and 
(2.23)  involves  an  inversion  of  the  matrix  shown  on  the  right  hand  side 
of  Equation  (2.23)  at  each  discrete  frequency  and  is  known  as  the  direct 
frequency  response  method  of  analysis.  The  inversion  of  this  usually 
large  and  complex  matrix  is  computationally  intensive  even  with  the  use 
of  conventional  elimination  procedures  [9]. 

An  attractive  alternative  to  using  a  direct  analysis  is  to  use  the 
normal  mode  approach  mentioned  earlier  on.  In  the  application  of  this 
method,  it  is  first  required  to  perform  an  eigenvalue  analysis  of  the 
structural  system  to  determine  its  natural  frequencies  and  mode  shapes. 
The  orthonormal i zed  mode  shapes  or  eigenvectors  are  used  to  form  a  modal 
transformation  matrix  R  which  is  used  to  uncouple  the  system  of  equations 
of  motion  into  normal  modes  or  generalized  coordinates.  The  main  advan¬ 
tage  gained  by  dealing  with  the  uncoupled  equations  is  that  the  complex 
frequency  response  function  matrix  is  now  a  diagonal  matrix  and  hence  the 
determination  of  its  inverse  for  use  in  the  relation  (2.22)  becomes  a 
trivial  computational  operation.  The  only  price  paid  is  that  the  modal 
response  obtained  must  be  transformed  back  to  the  original  global  degrees 
of  freedom  of  the  finite  element  model  of  the  structure.  A  restriction 
of  the  normal  mode  procedure  is  that  the  damping  matrix  must  be  such 
that  the  same  transformation  matrix  R  that  diagonalizes  the  mass  and 
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stiffness  matrices  will  do  the  same  for  the  damping  matrix.  For  this  to 
be  so,  the  damping  matrix  must  be  some  form  of  a  linear  combination  of 
the  mass  and  stiffness  matrices.  This  is  known  as  Rayleigh  or  propor¬ 
tional  damping  and  is  considered  to  be  an  adequate  representation  of  the 
damping  in  many  structural  applications.  It  is  also  conmon  practice  to 
assume  a  modal  damping  ratio  for  each  of  the  vibration  modes  used  for  the 
response  computations.  Another  disadvantage  of  the  modal  method  is  that 
only  some  of  the  modes  are  used  in  the  computation  of  the  responses  and 
in  some  cases  the  omitted  modes  could  contribute  significantly  to  the 
response  even  though  their  frequencies  lie  outside  the  range  of  the 
exciting  forces  [9],  This  problem  can  be  overcome  by  using  the  residual 
flexibility  method  which  is  not  implemented  in  this  work.  In  spite  of 
these  shortcomings,  the  normal  mode  approach  provides  an  elegant,  power¬ 
ful,  and  cost  effective  method  for  dealing  with  the  random  response  of 
discretized  structures.  As  such,  it  has  been  popularly  applied  in 
research  studies  and  commercial  finite  element  programs  [10,11,12]. 

Consider ,  once  again,  the  discretized  matrix  equation  of  motion 
given  in  (2.21).  The  natural  frequencies  (<or)  and  mode  shapes  (<j>r) 
of  the  system  are  first  computed  to  determine  [R]  using  M  and  K  which  are 
determined  by  the  finite  element  program.  Next,  the  assembled  cross- 
spectral  density  matrix  [Spp(to)]  must  be  determined  as  a  function  of 
the  forcing  frequencies  of  interest  from  the  user  specified  auto-  and 
cross-spectral  densities  of  the  applied  point  loads  and/or  distributed 
loads. 

Let  a  transformation  between  the  global  degrees  of  freedom  X  of  the 
finite  element  model  and  the  NM  modal  (general zed)  coordinates  q  be 
defined  such  that 

*  =  5  9  (2.22) 

where  R  is  the  mass-orthonormal ized  modal  matrix  described  earlier.  Sub¬ 
stituting  the  transformation  relation  (2.22)  in  (2.21)  and  pre-multiply- 
ing  both  sides  of  the  latter  by  the  transpose  (RT)  of  R  gives: 
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RTMRq  +  RTCRq  +  RTKRq  =  RT  F(t)  (2.23) 

in  which  the  nature  of  R  is  such  that 

RT  M  R  =  I  (2.24a) 

RT  C  R  =  [2£rwr]  (2.24b) 

RT  K  R  =  [u>2]  (2.24c) 

fr(t)  =  RTF(t)  (2.24d) 

R  =  [R]  =  [f4>x}  {♦,}  ...  {<JiNM}]  (2.24e) 


In  the  above,  I  is  the  (NMxNM)  identity  matrix,  [2^r«r  ]  and  [top  are 
(NMxNM)  diagonal  matrices  in  which  €p  is  the  equivalent  modal  damping  for 
the  rth  mode,  o>r  is  the  natural  frequency  corresponding  to  mode  r  out  of 
a  possible  total  number  of  NM  modes,  fp(t)  is  the  random  modal  force, 
and  {<t>r},  r=l,2,...,NM,  are  the  eigenvectors  of  the  NM  modes.  Thus, 
Equation  (2.23)  may  be  alternatively  written  as 


I  Qr  +  [2?ra)r]qr  +  [o>r2]qr  =  fp(t),  r=l,2,...NM.  (2.25) 

In  the  normal  mode  procedure  adopted  in  this  work,  we  shall  solve 
Equation  (2.25)  in  a  probabilistic  sense  to  obtain  the  random  modal 
displacement  response  and  convert  back  to  the  response  referred  to  the 
original  degrees  of  freedom  using  the  transformation  relation  defined  in 
(2.22).  This  is  because  the  uncoupled  set  of  Equations  (2.25)  is  easier 
to  handle  than  the  original  Equation  (2.21)  because  of  reasons  we  have 
previously  described. 
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The  first  step  is  to  determine  the  cross  spectral  density  matrix  of 
the  modal  force  [Sff(w)]  from  the  cross  spectral  density  matrix  of  the 
applied  force  referred  to  the  original  global  degrees  of  freedom  using 
the  relation: 


[Sff(w)]  -  [R]T  [SpF(w)]  [R]  (2.26) 
(NMxNM)  (NMxNS)  (NSxNS)  (NSxNM) 

This  relation  follows  directly  from  the  linear  transformation  relation 
between  f(t)  and  F(t)  defined  in  Equation  (2.24d).  An  important  property 
of  cross  spectral  density  matrices  of  real  processes  is  that  they  are 
Hermitian,  that  is. 


CSpF(w)]  -  [Spp(w)]T  (2.27a) 
or 

sJp(u)  =  SjjjL  (u)  (2.27b) 

The  Hermitian  property  of  complex  matrices  is  a  generalization  of  the 
symmetry  property  of  real  matrices.  Thus  the  cross-spectral  density  of 
the  input  excitation  to  the  modal  system  described  in  Equation  (2.25)  is 
now  known. 

The  second  step  is  to  compute  the  cross-spectral  density  of  the 
modal  displacement  response  in  Equation  (2.25)  using  the  complex  frequen¬ 
cy  response  functions  of  this  modal  system  and  the  input-output  relations 
expressed  in  (2.22).  Thus,  the  cross-spectral  density  of  the  random 
modal  displacement  response  is  given  by 


[Sqq(w)]  =  CH(ico)]  [Sff(w)]  [H(iw)]T* 


(2.28) 
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Equation  (2.28)  represents  the  frequency-domain  probabilistic  description 
of  the  solution  to  the  modal  system  of  Equations  (2.25).  We  emphasize, 
once  again,  that  the  computational  advantage  of  solving  Equation  (2.25) 
instead  of  (2.21)  lies  in  the  fact  that  the  complex  frequency  response 
function  in  (2.28)  is  a  diagonal  matrix  thus  making  it  easier  to  deal 
with.  The  frequency  response  function  [H(iu>)]  matrix  is  computed  accord¬ 
ing  to  the  type  of  damping  used. 

For  the  present  work,  the  damping  is  classified  into  three  cate¬ 
gories.  The  first  type  of  damping  requires  the  specification  of  the 
damping  ratio  £r,  for  each  mode  r  of  the  NM  modes  involved.  This  is  a 
ratio  of  viscous  damping  to  the  critical  damping.  Denoting,  the  element 
(r,r)  of  the  diagonal  matrix  [H(iw)]  by  Hrr(iw),  then  for  this  type  of 
modal  damping  for  any  forcing  frequency  w 

Hrr(i<*>)  =  {(«r2  "  w2)  +  (2^ro)r(u)i}-1  (2.29a) 

The  second  type  of  damping  is  the  proportional  damping  or  Rayleigh  damp¬ 
ing  described  earlier.  For  this  case,  the  damping  matrix  is  assumed  to 
be  of  the  form 

C  =  aMM  +  a^K  (2.29b) 

where  and  a<  are  real  constants. 

For  this  case,  Hrr(ico)  is  given  by 

Hrr(1u)  =  {(wr2-«2)  +  (Qm+VY2)  ui}-1  (2.29c) 

The  third  type  of  damping  option  is  structural  damping  for  which  the 
damping  force  is  represented  as  a  complex  stiffness  term  in  the  form: 


2.17 


CX  =  igKX 


(2.29d) 


where  g  is  the  structural  damping  factor. 

Unlike  modal  damping  and  Rayleigh  damping,  structural  damping  is  a 
non-vi scous  form  of  damping  that  is  usually  employed  for  modelling 
internal  energy  dissipation  per  cycle  (due  to  internal  friction)  for 
structures.  Thus,  for  this  third  case,  Hpp(ito)  is  given  by: 

HrrOw)  =  {(Wp2-w2)  +  (gu>r2)i}-1  (2.29e) 


All  these  three  damping  options  have  been  implemented  with  the  user 
defining  his/her  choice  by  means  of  the  parameter  IDAMP  which  may  assume 
values  of  1,  2,  or  3. 

The  final  step  is  the  transformation  of  the  cross-spectral  density 
of  the  modal  response  to  the  displacement  response  referred  to  the 
original  global  degrees  of  freedom,  using  the  relation: 

Csxx(w):  =  CR]  CSqq(«)3  [R]T  (2.30) 

(NSxNS)  (NSxNM)  (NMxNM)  (NMxNS) 

Again,  Equation  (2.30)  is  a  direct  consequence  of  the  linear  transforma¬ 
tion  relation  between  X  and  q  defined  by  the  expression  (2.22).  All  the 
spectral  density  matrices  in  (2.28)  and  (2.30)  satisfy  the  same  Hermitian 
property  described  in  (2.27).  Thus,  Equation  (2.30)  is  the  frequency- 
domain  probabilistic  description  of  the  random  displacement  response. 
Formally,  a  knowledge  of  the  cross-spectral  density  matrix  of  the 
displacements  completes  the  solution  process.  However,  in  practice,  a 
knowledge  of  other  statistical  properties  of  displacements  or  those  of 
other  related  processes  are  usually  of  interest.  These  related  processes 


include  derivative  processes  such  as  velocities  and  accelerations  and 
secondary  processes  such  as  stresses  and  strains  that  are  directly 
(linearly)  related  to  displacements. 

The  cross-spectral  densities  of  velocities  and  accelerations  can  be 
obtained  from  those  of  displacements  by  using  the  property  of  derived 
processes.  Denoting  the  cross  spectral  density  matrices  of  velocites  and 
accelerations  as  [Svv(w)]  and  [Saa(w)]  respectively,  then  the 
following  relations  are  used  for  their  computation: 

[Svv(u>)]  =  «2CSxx((o)]  (2.31) 

[SM(co)]  =  w“CSxx(a))]  (2.32) 

For  stresses  and  strains  (referred  to  as  secondary  responses  in  this 
work)  the  procedure  is  not  as  straightforward  as  those  given  in  equation 
(2.31)  and  (2.32)  above.  This  is  because  neither  the  stress  nor  the 
strain  is  a  time  derivative  of  the  displacement.  Although,  for  linear 
finite  element  analysis,  the  stresses  and  strains  are  related  to 
displacements  through  a  linear  transformation,  some  additional  effort  is 
required  for  the  determination  of  the  transformation  matrices.  Once 
these  transformation  matrices  are  determined,  the  rule  used  for  the 
computation  of  the  cross-spectral  densities  of  modal  force  from  those  of 
the  global  nodal  forces  can  then  be  applied. 

Let  [T]  be  the  transformation  matrix  relating  the  stress  {o}  at  a 
given  set  of  points  to  the  modal  displacements  {q}  and  let  [W]  be  the 
corresponding  matrix  for  the  strains  {e}.  Then 

{0}  =  [T]{q}  (2.33) 


and 


{£}  =  [W]{q] 


(2.34) 
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A  typical  component  T-j j  of  the  matrix  [T]  represents  the  stress  at 
point  1  due  to  a  normalized  oscillation  in  the  oscillation  mode  j.  m 
other  words,  therefore,  [T]  is  the  modal  stress  transformation  matrix 
corresponding  to  the  points  and  stress  components  of  interest  in  the 
structure.  The  same  goes  for  the  definition  of  the  components  of  [W] 
Comparing  Equations  (2.33)  and  (2.34)  with  the  transformation  relation 
(2.22),  1 1  is  easy  to  see  that  the  [T]  and  [W]  play  the  same  role  for 
stresses  and  strains  as  the  modal  matrix  [R]  does  for  displacements. 
Once  [T]  and  [W]  have  been  determined,  the  cross  spectral  density 
matrices  of  stresses  and  strains  are  obtained  from  the  relations: 

CS00(o»]  =  CT]  HSqq(u>)]  CI]T  (2.35) 

and 

CSee<")3  -  CSqqC»)3  [Wf  (2  36) 

WHh  Equations  (2.30),  (2.31),  (2.32),  (2.35),  and  (2.36),  the  cross 
spectral  densities  of  anj  of  the  random  response  quantities  that  may 
potential ly  be  of  Interest  In  a  random  vibration  analysis  may  be  comput¬ 
ed.  However,  In  practical  applications,  other  statistical  properties  of 
the  random  responses  may  be  additionally  required  or  indeed  preferred  to 
the  cross  spectral  densities.  In  what  follows,  a  brief  description  of 
these  other  properties  are  presented- 

Let  X1=Xl(t)  and  X2=X2(t)  be  two  random  processes  or  components  of  a 
random  vector  process.  The  covariance  of  Xx  and  X2  denoted  Cov(XxX2)  or 
cXxX2  is  defined  as 

Cov(XxX2)  =  CXiX2  =  E[[xx(t)  -  lij(t)].[X2(t)  -  y2(t)J]  (2.37) 

where  jjx  and  p2  are  the  mean  values  of  Xx  and  X2  respectively  and  the 
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operator  E  denotes  the  mathematical  expectation  defined  earlier  on.  For 
zero  values  of  ^  and  p2, 

Cx  x  =  ECXiX,]  (2.38a) 

and  for  the  special  case  where  Xx  and  X2  are  the  same  random  process 
(i.e.  X,(t)  =  X2(t)),  then 

Cv  y  -  C„  x  =  E[XxXJ  =  E[Xx 2]  (2.38b) 

AiA2  A i Ax 

Thus,  the  covariance  of  Xj  and  itself  usually  called  the  autocovariance 
(or  simply  variance)  of  Xt  is  nothing  but  the  mean  square  value  of  Xj 
provided  X!  is  a  zero-mean  process  which  we  assume  throughout  this  work. 
Let  {y}  be  a  random  vector  process  of  order  n.  Then  the  covariance  of 
{y}  is  an  nxn  matrix  defined  as 

[Cyy]  -  E[{yWy}T] 
or 

E[yi2]  ECyxy2]  ...  E[y1yn] 

E[y22] 

Sym. 

E[y=] 

It  can  be  seen  from  Equation  (2.39b)  that  for  zero-mean  processes,  the 
diagonal  terms  of  the  covariance  matrix  give  the  mean  square  responses. 

The  cross  and  auto  correlation  functions  are  related  to  the  cross 
and  auto  spectral  density  functions  through  the  well-known  Wiener- 
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Khintchine  relations  which  were  presented  earlier  (Equations  (2.7)  and 
(2.8)).  The  cross-correlation  function  between  x  and  y  is  given  by 


Rxy(x)  =  70  Sxy(w)eia,t  dw  (2.40) 

Thus,  with  a  knowledge  of  Sxy(w),  Rxy(t)  may  be  computed  by  taking 
the  inverse  discrete  Fourier  transform  (DFT)  of  Sxy(w)  using  the  FFT 
algorithm.  For  auto-correlation  functions,  in  particular 

ao  oo 

Rxx(t)  =  /  sxx(w)  e1U)Tdo)  =  2/  Sxx(w)  cos  «i  dw  (2.41a) 

-®  o 


or 


oo 

Rxx(t)  =  /  GXX(U))  cos  wx  dw 

_ oo 


(2.41b) 


where 

Gxx(o>)  =  2Sxx(w)  (2.41c) 

is  called  the  one-sided  power  spectral  density  which  has  zero  values  for 
negative  values  of  the  circular  frequency,  w. 

Recalling  that 

E[x!].«„Wl,.0-W°l  (2'423) 

it  is  easy  to  see  from  (2.41)  that 


E[x2]  =  /  S  (w)d«  =  2/  S  (w)  dw 

-CO  0 


(2.42b) 


Similarly,  since  E[x(t)y(t«)]  =  Rxy(t) ,  and 

Rxy(T) I t=0  =  £tx(t)y(fc)3  =  Cov(xy)  »  Cxy  (2.43a) 

it  follows  that 


(2.43b) 


Thus,  the  off-diagonal  terms  of  the  covariance  matrix  defined  in  (2.39b) 
may  be  computed  from  the  cross-spectral  densities  using  (2.43)  while  the 
diagonal  terms  may  be  computed  from  the  auto  spectral  densities  using 
(2.42).  The  diagonal  terms  popularly  referred  to  as  the  variances  or 
mean  square  responses  are  often  of  interest  while  off-diagonal  terms  may 
also  be  of  interest  in  reliability  studies. 

The  above  procedures  for  calculating  the  covariances  of  the 
responses  is  straightforward.  However,  it  is  time-consuming  from  a 
computational  implementation  point  of  view. 

Let  us  consider  a  random  vector  process  {y}  that  is  related  to 
another  random  vector  process  [z]  through  a  linear  transformation  matrix 
[T]  such  that 

(z)  =  CT][y}  (2.44) 
The  cross  spectral  density  matrices  of  [z]  and  £y)  are  related  through 

CSz(cu)]  =  [T][Sy((0)][T]T  (2.45) 


Similarly,  the  covariance  matrix  of  {z}  and  [y]  are  related  via  the 
expressions 
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[Cz]  =  [T][Cy][T]T  (2.46) 

Suppose  the  cross  spectral  density  [Sy(w)]  of  {y}  is  known  and  the 
covariance  of  (z)  is  to  be  computed.  Then  two  alternative  approaches  are 
available.  In  the  first,  [Sz(w)]  is  computed  for  all  the  applied  fre¬ 
quencies  using  the  relation  (2.45).  Then  [Cz]  is  determined  by  using 
the  matrix  form  of  the  relation  (2.43),  namely 


[C  T  =  2/  [S  (w)]  du>  (2.47) 

z  0  z 

In  the  second  approach,  the  covariance  matrix  [Cy]  of  (y)  is  computed 
using  the  relation 


[Cv]  =  2/  [S  (o))]  du>  (2.48) 

y  q  y 

and  then  the  covariance  matrix  of  {z}  is  computed  using  equation  (2.46). 

It  can  be  seen  that  the  first  approach  involves  the  calculation  of  the 
triple  matrix  product  on  the  right  hand  side  of  (2.45)  for  each  of  the 
applied  forcing  frequencies.  For  a  total  number  of  NAF  applied  frequen¬ 
cies,  therefore,  the  triple  matrix  product  has  to  be  evaluated  NAF  times 
and  so  the  procedure  is  computationally  intensive.  In  the  second 
approach,  however,  the  integration  in  Equation  (2.48)  is  first  performed 
and  then  the  triple  matrix  product  in  (2.46)  is  evaluated  only  once. 
Thus,  the  first  approach  involves  one  numerical  integration  and  NAF 
evaluations  of  a  triple  matrix  product  while  the  second  involves  one 
numerical  integration  and  the  evaluation  of  only  one  triple  matrix 
product.  The  second  approach  is  much  more  preferred,  especially  when  the 
structural  model  size  is  large  or  when  there  are  more  than  a  few  forcing 
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frequencies  because  of  the  significant  savings  in  computation  time.  The 
only  situation  in  which  the  first  procedure  would  be  preferable,  of 
course,  is  when  the  cross  spectral  density  [Sz(w)]  has  been  previously 
computed  or  is  otherwise  known. 

In  terms  of  the  response  quantities  of  interest  in  this  work,  the 
above  discussion  of  the  computation  of  covariances  of  responses  trans¬ 
lates  to  the  following.  For  displacements,  for  example,  if  the  spectral 
densities  are  not  of  interest  but  say  the  mean  squares  are,  then  they  are 
best  computed  by  first  computing  the  covariance  matrix  of  modal  displace¬ 
ments  from  the  corresponding  cross  spectral  density  matrix,  that  is 


[Cqq]  =  2/o[Sqq(o»]  do)  (2.49) 

It  will  be  recalled  that  [Sqq(u))]  has  been  previously  computed  as  the 
second  step  in  the  random  response  analysis,  Equation  (2.28).  Then  the 
covariance  matrix  of  displacements  [Cxx]  is  computed  from 

[Cxx]  =  [R][CQq][R]T  (2.50) 

The  diagonal  terms  of  [Cxx]  represent  the  mean  square  displacements. 
Similar  procedures  apply  to  other  response  quantities.  For  example,  for 
stresses,  the  covariance  matrix  is  computed  from  the  relation 

[c0<J:  -  :TKcqq:[T:T  <2.5» 

and  the  diagonal  terms  of  [Coa]  give  the  mean  square  stresses.  Higher 
order  spectral  moments  of  any  of  the  random  response  quantities  may  also 
be  computed  by  using  the  cost  effective  second  approach.  This  was  used, 
for  example,  by  Karadeniz  et  al  [13]  for  the  spectral  analysis  of  the 
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probabilistic  reliability  analysis  of  the  fatigue  limit  state  of  gravity 
and  jacket-type  structures. 

The  mean  square  response  is  a  very  important  statistical  property 
because  it  can  be  used  for  calculating  the  probability  of  failure  (and 
hence  the  reliability)  of  a  structural  system  once  the  probability 
density  function  of  the  random  response  process  is  known.  It  is  an 
interesting  fact  of  life  that  many  naturally  occurring  random  vibrations 
have  the  well-known  "bel 1 -shaped"  probability  distribution  function  known 
as  the  normal  or  Gaussian  probability  distribution.  Accordingly,  a 
random  process  that  possesses  this  type  of  distribution  is  known  as  a 
Gaussian  process.  An  important  property  of  linear  time-invariant  systems 
excited  by  Gaussian  random  process  is  that  the  resulting  response  is  also 
Gaussian  and  hence  possesses  the  same  probability  density  function  as  the 
excitation.  The  probability  density  function  of  the  random  process  x  has 
the  form 

p(x)  =  — - —  exp  [-(x-m)2/2o2]  (2.52) 

/2i  o 

where  m  is  the  mean  E[x]  and  o2  is  the  variance  E[(x-E[x]) *].  The  normal 
distribution  function  is  completely  characterized  by  the  two  parameters  m 
and  a2.  For  a  zero-mean  process  (i.e.  for  m=0) ,  o2  is  the  same  as  the 
mean  square.  Thus  the  mean  square  value  of  x  can  be  used  to  compute  the 
probability  that  x=x(t)  will  lie  in  a  certain  interval  (say  x^t)  and 
xa(t))  using  the  relation 


x2 

Prob[X!  <  x(t)  <  x2]  =  P(x2)  -  P(xx)  =  J  p(x)dx  (2.53) 

Xi 


Note  that  in  equation  (2.53),  P(x-j)  denotes  the  probability  that  x=x(t) 
will  not  exceed  the  value  x^ ,  that  is. 
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Prob[x(t)  <  x.]  =  /  p(x)dx 


(2.54) 


Thus,  if,  for  example,  failure  is  assumed  to  occur  if  x=x(t)  goes  outside 
the  range  x!<x(t)<x2,  then  from  (2.53),  the  probability  of  failure  is 
1  -  Prob  [Xl<x(t)<x2].  Consequently,  the  reliability  is  1  -  (probability 
of  failure)  =  Prob[x!<x(t)<x2].  If  failure  is  assumed  to  occur  if  x=x(t) 
exceeds  a  critical  value  x^,  say,  then  the  probability  of  failure  is 
Prob[x(t)>x.j]  =  1  -  Prob[x(t)<x.j]  which  is  readily  obtained  from  (2.54). 
In  this  case,  the  reliability  is  Prob[x(t)<x because  the  reliability, 
r,  is  always  related  to  the  probability  of  failure,  pf,  through 


r  =  1 


(2.55) 


For  the  cases  where  there  are  two  or  more  random  variables  or  for 
vector-valued  random  processes  such  as  we  have  in  the  analysis  of  multi- 
degree-of-freedom  systems,  the  first-order  probability  density  function 
is  no  longer  sufficient  to  describe  the  probability  structure  of  the 
Gaussian  process.  A  full  definition  of  the  multi-dimensional  Gaussian 
distribution  must  include,  not  only  the  first-order  probability  density 
function  given  by  (2.52),  but  also  corresponding  second  and  higher  order 
joint  densities.  For  the  case  of  two  jointly  Gaussian  random  variables  x 
and  y,  for  example,  the  second-order  probability  density  function  is 
[5]: 


P(x,y)  = 


2™xy(i-p2Xy) 


(x-mx)2  (y-mY)2  2pxv(x-mx)(y-my) 

Xy)  o2x  02y  OxOy 


(2.56) 


where  m  and  m  are  the  mean  values  of  x,  and  y,  o2  and  o*  are  the  var- 
x  y  x  y 

iances  of  x  and  y  and  p  is  a  correlation  coefficient  defined  as 

xy 


1 
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E[(x-m  ) (y-m  )] 


o  o 

x  y 


(2.57) 


and  called  the  normalized  covariance.  Notice  that  if  Px  =  0.  x  and  y  are 
statistically  independent  because  then  (2.56)  may  be  factorized  to  give 

P(x,y)  =  (—3 —  e-(x-mx2/2ox2)i,r — ! —  e-(y-my2/2av2) l 

^  V2io  y  Y  J 

x  y 


=  P(x)*p(y) 


(2.58) 


The  numerator  of  the  right  hand  side  of  Equation  (2.57),  it  will  be 
recalled,  is  a  typical  term  of  the  covariance  matrix  of  x  and  y.  it 
represents  an  off-diagonal  term  for  x*y  and  a  diagonal  term  for  x=y. 
Thus  it  can  be  seen  that  the  off-diagonal  elements  of  the  covariance 
matrix  are  also  important  In  the  application  of  higher  order  joint 
densities  for  the  computation  of  failure  probabilities  or  in  reliability 
analysis.  Zero-valued  off-diagonal  terms  of  the  covariance  matrix 
correspond  to  zero  correlation  and  the  higher-order  joint  density 
functions  become  simply  the  products  of  the  first-order  density  functions 
as  in  (2.58).  However,  for  MDF  system  excited  by  random  excitations  it 
is  not  uncommon  for  the  off-diagonal  terms  to  take  on  non-zero  values 
because  the  response  process  is  actually  a  vector  random  process  in  which 
the  components  are  like  correlated  random  variables.  This  is  the  motiva¬ 
tion  for  including  the  capability  for  the  computation  of  cross  spectral 
densities  and  cross  terms  of  the  covariances  of  responses  in  this  work. 
It  is  hoped  that  future  work  involving  the  application  of  advanced  reli¬ 
ability  analysis  techniques  will  benefit  intensely  from  the  broad  frame¬ 
work  that  has  been  established  during  the  course  of  this  first  phase  of 
the  contract. 


Another  statistical  property  of  response  that  Is  usually  of  interest 
in  probability  structural  dynamics  is  the  apparent  frequency  or  the 
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expected  (or  effective)  equivalent  frequency.  This  is  defined  as  the 
expected  value  of  the  number  of  zero  crossings  with  positive  slope  per 
unit  time  denoted  N0  or  E[N+(o)].  For  a  random  process  x-x(t)  whose 

corresponding  time-derivative  process  is  x=x(t),  it  is  mathematically 
defined  as: 


N0  ■  E[N+(o)]  .  | 


(2.59) 


where  ox  is  the  variance  of  x  and  ox  is  the  variance  of  the  corresponding 
rate  process  [14]].  For  a  zero  mean  process  then 


J  <*>2  S  (»)  du> 

N  =  i-  X 

o  2  it  - - 

/  sxx(w)  du> 


(2.60) 


In  the  relation  (2.60),  Sxx«o)  is  the  auto  spectral  density  of  x  and  so 
it  1S  easy  to  see  that  the  denominator  is  the  mean  square  of  x  while  the 
numerator  Is  the  mean  square  of  X.  In  view  of  the  discussions  previously 
given  concerning  the  efficient  computation  of  mean  square  responses  a 
knowledge  of  the  spectral  density  of  x  is  not  essential  in  the  computa¬ 
tion  of  its  apparent  frequency.  Both  the  denominator  and  the  numerator 
in  (2.60)  can  be  computed  directly  from  the  covariances  of  the  modal 
displacements  and  modal  velocities  respectively.  The  apparent  frequency 
is  very  useful  in  fatigue  analysis  and  so  the  capability  for  computing 

this  statistical  property  for  any  of  the  possible  physical  processes  has 
been  provided  in  this  work. 


3.1 


3.0  RANDOM  LOAD  DESCRIPTION 

The  analysis  procedure  presented  in  the  preceding  chapter  is  based 
on  the  assumption  that  the  cross  spectral  density  matrix  of  the  applied 
forces  referred  to  the  original  global  degrees  of  freedom  of  the  finite 
element  model  is  available.  However,  this  is  not  the  case  in  practice. 
Although  the  user  defines  the  spectra  and/or  cross  spectra  of  the  point 
loads  and/or  distributed  (i.e.  pressure)  loads  acting  on  the  structure,  a 
procedure  is  required  for  calculating  the  consistent  assembled  cross 
spectral  density  matrix  of  the  system  due  to  all  the  applied  random  loads 
in  a  finite  element  setting.  This  operation  is  analogous  to  the  computa¬ 
tion  of  the  assembled  consistent  load  vector  in  a  conventional  determin¬ 
istic  analysis  using  the  finite  element  method. 

The  primary  goal  of  the  present  work  is  the  development  of  the  capa¬ 
bility  for  the  random  response  analysis  of  naval  structures  travelling  in 
rough  seas  so  that  more  accurate  and  realistic  predictions  of  the  integ¬ 
rity  and  safety  of  such  structures  can  be  made  by  the  scientific  author¬ 
ity.  Accordingly,  the  main  source  of  loading  that  must  be  dealt  with  is 
the  random  pressure  loading  imparted  onto  the  hull  structure  due  to  the 
action  of  the  random  ocean  waves.  This  form  of  loading  represents  the 
most  significant  time-varying  loads  experienced  by  marine  structures 
[15].  Since  the  adequate  description  of  the  loading  is  one  of  the 
fundamental  steps  involved  in  the  analysis  of  any  engineering  structure, 
the  development  of  a  procedure  to  appropriately  describe  the  loading  in 
the  framework  of  a  probabilistic  finite  element  methodology  is  a  very 
important  component  of  the  work  presented  in  this  report. 

3*1  Wave  Loads  on  Marine  Structures 

The  computation  of  hydrodynamic  wave  forces  on  marine  structures  is 
a  difficult  task  because  it  involves  the  complexity  of  the  interaction  of 
waves  with  the  structure.  The  two  approaches  commonly  used  may  be 
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broadly  classified  into  deterministic  and  statistical.  Although  the 
statistical  approach  is  the  one  of  interest  in  this  work,  it  is  quite 
appropriate  to  give  a  brief  overview  of  the  essential  ingredients  of  each 
of  the  two  approaches  because  the  statistical  approach  is  not  completely 
independent  of  the  deterministic  one.  This  interdependence  will  become 
evident  as  we  progress  in  this  report. 

In  a  deterministic  approach,  an  appropriate  wave  theory  for  describ¬ 
ing  the  water  particle  kinematics  for  the  given  wave  condition  is  select¬ 
ed  and  the  hydrodynamic  forces  are  determined.  The  three  different 
methods  for  this  class  of  approach  include  the  Morison  equation,  the 
Froude-Kryl ov  theory,  and  the  diffraction  (or  potential  flow)  theory 
[16].  The  Morison  equation  assumes  the  force  to  be  a  linear  combination 
of  inertia  and  drag  forces.  The  inertia  and  drag  components  of  the  force 
involve  coefficients  which  are  determined  experimentally.  The  Morison 
equation  is  very  suitable  for  describing  the  hydrodynamic  loading  on 
slender  structures  since  the  drag  force  is  significant  for  such  applica¬ 
tions.  For  example,  the  in-line  force  per  unit  length  (p)  on  a  vertical 
slender  cylinder  in  unsteady  flow  may  be  represented  in  the  form  [16]: 

P  ■  cdp  \  unl“J  +  cra  p’  IT‘n  f3'11 

where  p  is  the  water  density,  D  is  the  diameter,  un  and  0n  are  the  water 
particle  velocity  and  acceleration,  respectively,  and  Cj  and  Cm  are 
the  drag  and  inertia  coefficients,  respectively. 

For  situations  in  which  the  drag  force  is  small  and  the  inertia 
force  predominates,  but  the  structure  is  still  relatively  small  compared 
to  the  water  wave  length,  the  Froude-Krylov  theory  may  be  applied.  This 
method  utilizes  the  incident  wave  pressure  and  the  pressure-area  method 
on  the  surface  of  the  structure  to  compute  the  force.  Because  the 
calculation  of  the  force  on  the  structure  is  performed  assuming  that  the 
structure  is  not  there  as  far  as  the  waves  are  concerned,  this  approach 
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has  limited  practical  applications.  The  merit  of  the  method  lies  in  its 
simplicity  and  the  possibility  of  closed-form  force  solutions  for  some 
symmetric  structures. 

For  large  structures,  the  characteristic  dimension  of  the  structure 
becomes  comparable  to  the  wave  length  and  so  the  presence  of  the 
structure  is  expected  to  alter  the  wave  field  in  the  vicinity  of  the 
structure.  If  the  structure  spans  a  significant  portion  of  a  wave 
length,  the  incident  waves  upon  arriving  at  the  structure  undergo 
significant  scattering  or  diffraction.  Therefore,  the  diffraction  of  the 
waves  from  the  surface  of  the  structure  should  be  taken  into  account  in 
the  evaluation  of  the  wave  forces.  This  method  is  generally  known  as  the 
diffraction  theory.  The  solution  involves  the  numerical  solution  of  the 
Laplace  equation  together  with  the  associated  boundary  conditions, 
namely: 


V2$ 


a2$  a2$  ,  a2$ 

■  8x2  3y2  9z2  - 


(3.2) 


Equation  (3.2)  assumes  the  basic  flow  to  be  oscillatory,  incompressible, 
and  irrotational  so  that  the  fluid  velocity  may  be  represented  as  the 
gradient  of  a  scalar  potential,  $  =  fi(x,y,z,t).  Details  of  this  and 
other  theories  are  very  well  elaborated  upon  in  the  monograph  by 
Chakrabarti  [16]. 


The  second  approach  to  the  characterization  of  loads  induced  by 
ocean  waves  is  the  utilization  of  statistical  description  techniques. 
This  approach  is  considered  more  realistic  in  view  of  the  fact  that  real 
sea  surface  profiles  are,  more  often  than  not,  highly  irregular  and  con¬ 
fused  in  nature,  and  nonrepeatable  in  time  and  space.  In  this  approach, 
the  ocean  waves  and  hence  the  associated  sea  surface  profiles  are  consid¬ 
ered  to  be  random.  The  surface  elevation,  n»  is  considered  to  be  a 
random  process  which  is  described  by  either  its  probabilistic  properties 
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(such  as  the  cummulati ve  distribution  function,  P(n),  and  the  probability 
density  function,  p(n)),  or  the  power  spectral  density  S^n(k,w)  usually 
referred  to  simply  as  the  spectrum.  Unlike  the  probabilistic-function 
description,  the  spectral  description  method  preserves  the  frequency 
content  of  the  wave  field  or  surface  elevation  profile  and  the  relative 
distribution  of  energies  at  different  frequencies.  Furthermore,  it  is 
naturally  compatible  with  the  frequency  domain  method  for  response 
analysis  which  was  discussed  extensively  in  Chapter  2.  As  such,  it  has 
found  wide  applications  iri  engineering  analysis.  Thus,  the  wavy  sea 
surface  profile,  q,  is  in  general  modelled  as  a  stochastic  field,  which 
for  simplicity,  is  usually  considered  to  be  stationary  in  time  and 
homogeneous  in  space. 

3.2  Ocean  Wave  Spectra 

The  most  appropriate  procedure  for  obtaining  a  spectrum  is  by 
measurement  at  the  site  under  consideration.  Unfortunately,  such  a 
spectrum  is  seldom  available,  especially  when  directional  spectra  are 
required.  As  an  alternative,  the  usual  procedure  is  to  choose  one  of  the 
theoretical  spectrum  models  available.  This  theoretical  wave  spectrum 
describes  short  term  wave  conditions,  since  spectrum  models  depend  on 
fetch,  wind,  and  other  meteorological  conditions  of  the  site. 

Several  mathematical  spectrum  models  are  available.  These  models 
are  generally  based  on  one  or  more  parameters,  for  example,  significant 
wave  height  and  wave  period.  A  few  of  these  spectrum  models  that  have 
been  applied  to  marine  structures  include  the  Pierson-Moskowltz  spectrum, 
the  ISSC  spectrum,  the  ITTC  spectrum,  and  the  JONSWAP  spectrum.  A  brief 
description  of  each  of  the  above  spectra  is  given  below. 

The  Pierson-Moskowitz  spectrum  (usually  referred  to  as  the  P-M 
model)  is  a  single-parameter  spectrum  based  on  the  significant  wave 
height  or  wind  speed.  It  was  developed  in  1964  from  the  results  of 
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analysis  of  wave  spectra  in  the  North  Atlantic.  It  has  since  been 
extensively  used  by  ocean  engineers  as  one  of  the  most  representative  for 
waters  all  over  the  world  [16,17].  The  P-M  spectral  model  describes  a 
fully-developed  sea  determined  by  the  wind  speed.  The  fetch  and  duration 
are  considered  infinite.  The  model  is  applicable  to  wave  records 
obtained  after  the  wind  has  blown  over  a  large  area  at  a  nearly  constant 
speed  and  direction  for  many  hours  prior  to  the  time  of  taking  the  wave 
record.  It  is  mathematically  represented  as: 

S(o>)  =  ag2ur5  exp  [-0.74  (^)"4]  (3.3) 

where  a=0.0081,  g  is  the  acceleration  due  to  gravity,  and  Uw  is  the 
wind  speed.  Equation  (3.3)  may  also  be  written  in  terms  of  the  frequency 
of  the  spectral  peak,  <o0 ,  namely: 

S(»)  =  a92w-5  exp  [-1.25  (j^)-4]  (3.4) 

Karadeniz  et  al  [13],  used  a  multi-directional  P-M  spectrum  to 
characterize  the  sea  state  in  their  analysis  of  gravity  and  jacket-type 
structures.  Their  model  has  the  form: 

Snn^w’^  =  l-TT  cos2($-$s)][ag2w_B  exp{4ag2Hs-2ur4}]  (3.5) 

In  equation  (3.5),  $s  is  the  principal  direction  of  the  sea  state  in 
which  individual  waves  with  direction  $  simultaneously  occur,  o>  is  the 
frequency  of  wave  motion,  Hs  is  the  significant  wave  height,  g  is  the 
acceleration  due  to  gravity,  and  a  is  a  constant  (0.008). 

In  1964,  the  International  Ship  Structures  Congress  (ISSC)  suggested 
a  slight  modification  to  the  Bretschneider  spectrum  [16]  in  the  form: 
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S(w)  =  0.1107  Hs2  *£  exp[-0.4427 (jjj  )4] 
in  which  the  peak  frequency,  w0,  is  related  to  w  through: 


(3.6) 


0)  =  1.296io0 


(3.7) 


The  International  Towing  Tank  Conference  (ITTC)  proposed  (in  1966, 
1969,  1972)  a  modification  of  the  P-M  spectrum  in  terms  of  the 
significant  wave  height  and  zero  crossing  frequency,  wz.  The 
mathematical  representation  of  the  ITTC  spectrum  is  in  the  form: 


S(<d)  =  ag2or5 


exp[-«^ 


(3.8) 


where  a  = 


0.0081 

k4 


(3.9) 


and  k  - 

3.54u>z 

in  which  a  =  i/m^  =  4 


(3.10) 


(3.11) 


In  the  above,  0  is  the  standard  deviation  (r.m.s.  value)  of  the  water 
surface  elevation,  mn  is  the  nth  moment  of  the  spectrum  defined  as 


(3.12) 


and  the  average  zero  crossing  frequency  or  apparent  frequency  (defined  in 
Chapter  2),  wz  is: 
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The  JONSWAP  spectrum  was  developed  during  a  joint  North  Sea  Wave 
Project  from  which  its  name  is  derived.  The  formula  for  this  spectrum  is 
given  in  [16]  as: 


S(«) 


ag2w-5 


expC-1.25  <%-)-* 


(3.14) 


in  which  y  =  peakedness  parameter,  t  =  shape  parameter  (t  for  (o<w0  and 
Tb  for  w>u)0).  in  terms  of  the  significant  wave  height,  Hs,  and  the 

peak  frequency,  w0,  an  approximate  expression  for  the  JONSWAP  spectrum 
is: 


S(«) 


* 

a  H 


^  exp[-1.25(^)-]  Y 


exp 


where 


(3.15) 


*  _  _ 0.0624 _ 

0.230  +  0.0336y  -  0.185  (1.9+y)"1  (3.16) 

For  Y-l,  the  value  for  a*  is  0.312,  and  Equation  (3.15)  reduces  to  the 
P-M  spectrum.  The  JONSWAP  spectrum  applies  to  wind  generated  waves  under 
conditions  of  limited  fetch  and  homogeneous  wind  fields. 

A  host  of  other  spectrum  models  are  available  and  Chakrabarti  [16] 
has  a  quite  elaborate  description  of  these  models.  The  choice  of  a 
spectrum  model  in  the  design  or  analysis  of  marine  structures  is  up  to 
the  designer  or  analyst  and  depends  on  the  particular  application. 

All  the  spectrum  models  described  above  are  unidirectional  frequency 
spectra  which  do  not  account  for  the  directional  distribution  of  wave 
energy.  Since  ocean  waves  are  multi-directional  in  practice,  a  more 
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accurate  description  of  random  seas,  and  hence  sea  surface  elevation 
profiles,  requires  information  on  the  directional  dispersion  of  wave 
energy.  In  other  words,  the  sea  surface  elevation  is  a  random  function 
of  both  space  and  time  thereby  making  it  necessary  to  describe  the  power 
spectrum  in  terms  of  both  frequency  (w)  and  wave  number  (k).  Measure¬ 
ments  of  directional  wave  spectra  (or  wave  number  spectra)  are  consider¬ 
ably  more  difficult  than  those  of  frequency  spectra.  While  frequency 
spectra  can  be  obtained  from  a  wave  record  at  a  single  point,  directional 
wave  spectra  require  simultaneous  recording  of  several  wave  components 
[18].  A  few  measuring  techniques  in  respect  of  the  latter  case  are 
discussed  in  [18],  but  sources  of  measured  directional  spectra  are 
considered  to  be  rare  [17]. 

Some  methods  for  obtaining  directional  wave  spectra  discussed  by 
Ochi  [19]  include: 

(i)  obtaining  the  sea  surface  elevation  over  an  area  by  means  of 
stereophotographs; 

(ii)  analyzing  data  obtained  from  one-dimensional  arrays,  two- 
dimensional  arrays,  and  multielement  arrays; 

(iii)  obtaining  the  directional  information  from  wave  elevations, 
slope,  and  curvature,  as  recorded  by  a  floating  buoy. 

Nwogu  [20]  also  presented  a  procedure  for  estimating  directional  wave 
spectra  from  an  array  of  wave  probes  based  on  the  Maximum  Entropy  Method 
(MEM). 

Information  on  wave  directionality  is  especially  important  for 
predicting  wave-induced  forces  on  ships  and  floating  structures  because 
the  forces  are  associated  with  coupled  motions  induced  by  waves  travel¬ 
ling  in  various  directions  [19].  However,  because  of  the  difficulty  in 
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measuring  directional  spectra,  it  is  usual  to  evaluate  directional  wave 
spectra  by  applying  a  formulation  representing  the  energy  spreading 
function  to  a  point  (i.e.  unidirectional)  spectrum.  This  is  considered 
to  be  particularly  convenient  for  evaluating  the  responses  of  ships  and 
marine  structures  in  a  seaway  [19],  for  which  correlations  among  the  six 
degrees  of  freedom  have  to  be  considered. 

According  to  Nwogu  [20],  the  concept  of  a  directional  spreading 
function  to  characterize  a  wave  field  should  only  be  applicable  to  a 
spatially  homogeneous  wave  field  with  no  correlation  of  the  wave  compon¬ 
ents  travelling  in  different  directions.  In  [19],  the  MEM  was  reported 
to  have  yielded  a  spreading  function  of  each  frequency  band  which  is 

consistent  with  cross-spectra  of  the  measured  water  surface  elevation 
time  series. 


A  directional  wave  spectrum  is  usually  represented  in  the  form: 
S(w,0)  =  S(w)D(u,0)  {3>17) 

in  which  D(u>,0)  is  a  nonnegative  angular  spreading  function  describing 
the  distribution  of  wave  energy  over  direction,  and  has  to  satisfy: 


/^D(a),0)d0  =  1 


(3.18) 


In  Equation  (3.17),  the  wave  frequency  (w) 
(k)  by  the  linear  dispersion  relation, 
unidirectional  frequency  spectrum. 


is  related  to  the  wave  number 
and  S(w)  is  the  conventional 


Available  idealized  spreading  function  formulations  include  [19]  the 
cosine-square  formulation,  the  SWOP  formulation,  the  Mitsuyoisu  formula¬ 
tion,  and  the  Hasselmann  formulation.  The  SWOP  energy  spreading  func¬ 
tion,  for  example,  was  derived  from  analysis  of  data  obtained  during  the 
Stereo  Wave  Observation  Project  (SWOP),  and  has  the  form: 
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!—  (1  +  acos29  +  bcos49)  for  —?  <0<? 

77  22  (3.19) 

0  otherwise 

where 

a  =  0.50  +  0.82  exp  (-  |  a>4)  (3.20a) 

b  =  0.32  exp  (-  |  io4)  (3.20b) 

w  =  dimensionless  frequency  =  ~  (3.20c) 

and  U  is  the  wind  speed  at  19.5  m  above  the  sea  surface.  It  may  be  noted 

that  the  multi-directional  version  of  the  P-M  spectrum  used  by  Karadeniz 

et  aT.  [13],  that  is  Equation  [3.5],  is  a  cosine-square  formulation  of 
the  spreading  function.  This  is  extremely  simple  to  use  but  it  is  a 
function  of  neither  frequency  nor  wind  speed. 

Recently,  Juszko  [21]  employed  various  directional  parameterizations 
for  modelling  high-resolution  directional  wave  spectra  under  contract  to 
the  scientific  authority.  These  wave  spectra  were  obtained  from  a  WAVEC 
buoy  moored  on  the  Grand  Banks  in  1984.  The  ten-parameter  directional 
spectum  presented  has  the  form: 


Snn(<d,0)  =  4  1  (4^1±1  Wm-j4)^7^2  e 

4  1=1  4  1  _ 


COS 


2Pi(e-emi) 


r(*i)» 


4Ai+l 


(3.21) 


The  approach  taken  in  the  work  of  Reference  [21]  was  to  extend  the  Ochi 
and  Hubble  [22]  six-parameter  amplitude  model  to  include  various  repre¬ 
sentations  of  the  directional  component  through  the  use  of  a  mean  direc¬ 
tion  and  a  directional  spread  parameter.  The  model  was  reported  to  be 
capable  of  reproducing  over  90%  of  the  data  records  and  so  will  be  the 
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top  candidate  for  consideration  as  a  suitable  choice  of  the  sea  surface 
profile  in  this  work.  Further  discussions  on  directional  wave  spectra 
can  be  found  in  the  monograph  by  Price  and  Bishop  [23]. 

3.3  Determination  of  Random  Pressure  Loads 

An  important  step  regarding  the  probabilistic  descriptions  of  the 
loading  of  marine  structures  is  the  selection  of  an  appropriate  power 
spectrum  model  for  the  ocean  waves  as  discussed  in  the  preceding  sec¬ 
tion.  The  next  critical  step  is  the  utilization  of  this  information  for 
the  determination  of  the  corresponding  random  pressure  loads  experienced 
by  such  structures.  The  general  procedure  is  to  consider  the  spectrum  of 
the  sea  waves  or  sea  surface  elevation  as  the  input  process  and  seek  to 
determine  the  corresponding  force  cross  or  auto  spectra  densities  output 
with  the  aid  of  the  classical  relations: 

2n  * 

Spp(w)  =  J  Hpn  (1«.0)Hpn(i«,0)Snn(w,e)  d0  (3.22) 

in  which  SPP  is  the  spectral  density  of  the  pressure,  0  represents  the 
direction  of  the  individual  waves  making  up  the  sea  state,  Hpn 
represents  the  transfer  function  (or  frequency  response  function)  from 
sea  surface  elevation  to  wave  pressure,  and  Hj$n  denotes  the  complex  con¬ 
jugate  of  Hp^.  This  approach  has  been  used  by  Karadeniz  et  al  [13], 
for  example,  in  the  spectral  fatigue  analysis  of  offshore  platforms.  In 
their  work,  the  linear  wave  theory  and  Mori  son's  formula  were  applied  to 
the  undisturbed  wave  to  determine  these  transfer  functions.  They  also 
applied  a  correction  based  on  diffraction  theory  for  the  higher  wave 
frequencies.  It  is  in  the  determination  of  these  transfer  functions, 
either  experimentally  or  computationally,  that  the  classical  theories  of 
deterministic  hydrodynamic  wave  analysis  become  important  and  relevant  in 
a  probabilistic  approach. 


The  determination  of  the  transfer  functions  Hprjfiw.e)  whose 
magnitude,  in  the  terminology  of  hydrodynamcs,  is  more  popularly  referred 
to  as  the  response  amplitude  operator  (RAO)  is  needed  in  the  application 
of  Equation  (3.22)  which  is  valid  for  all  linear  or  linearized  systems. 
For  pressure  loads  which  are  of  interest  here,  the  RAO  is  formally 
defined  as  the  response  amplitude  per  unit  wave  height.  However,  it  is 
usually  more  convenient  in  hydrodynamics  to  describe  the  RAO  as  the 
response  amplitude  per  unit  wave  amplitude.  In  the  computation  of  RAO, 
the  waves  are  considered  regular  and  a  sufficient  number  of  frequencies 
(and  directions  in  case  of  directional  wave  spectra)  are  chosen  to  cover 
the  entire  range  of  frequencies  (and  directions  or  wave  headings)  covered 
by  the  ocean  wave  spectrum. 

In  the  current  work,  a.  program  called  PRECAL  is  used  for  the  compu¬ 
tation  of  the  transfer  functions.  This  program  is  one  of  the  three 
HPCFEM  (acronym  for  Hydrodynamic  Pressure  Calculation  for  a  Finite 
Element  Model)  suite  of  programs  that  is  available  to  the  scientific 
authority  as  a  member  organization  of  the  NSMB  Co-Operative  Research. 
The  purpose  of  this  suite  of  programs  is  the  calculation  of  the 
hydrodynamic  pressure  on  the  wetted  surface  of  a  ship  moving  in  regular 
sinusoidal  waves  for  input  into  a  finite  element  model.  The  first 
program  in  the  suite  called  HYDMES  is  used  for  the  generation  of  the 
hydrodynamic  mesh  data  from  a  given  ship  geometry  database.  The  second 
program  PRECAL  then  uses  the  hydrodynamic  mesh  data  to  calculate  the 
hydrodynamic  pressure  acting  on  the  wetted  surface  of  the  ship  caused  by: 

i.  Oscillatory  motion  in  regular  sinusoidal  waves  for  infinite 
depth; 

ii.  Forward  motion  in  regular  sinusoidal  waves  for  infinite  depth. 

The  hydrodynamic  pressure  at  a  point  on  the  wetted  hull  surface  per  unit 
wave  amplitude  is  derived  based  on  a  three-dimensional  linear  potential 
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theory  in  which  the  effects  of  the  incident  wave,  diffracted  wave,  radi¬ 
ated  wave,  and  the  resultant  vertical  displacement  of  that  point  are  all 
accounted  for.  The  fluid  is  assumed  to  be  inviscid  and  incompressible 
while  the  flow  is  taken  to  be  irrotational .  The  third  program  in  the 
suite,  FINMES,  is  the  Finite  Element  Mesh  program  which  provides  the  user 
with  the  hydrodynamic  pressure  values,  interpolated  from  the  values 
calculated  in  PRECAL,  for  either  a  finite  element  model  or  for  user 
specified  points.  It  must  be  pointed  out  that  the  analytical  basis  of 
all  computations  performed  in  the  entire  HPCFEM  suite  of  programs  is 
completely  deterministic  and  as  such  has  no  element  of  randomness  whatso¬ 
ever.  Further  details  of  the  theoretical  basis  and  the  operation  of  the 
HPCFEM  suite  of  programs  can  be  found  in  the  publications  of  the  Research 
and  Development  Division  of  the  American  Bureau  of  Shipping  and  allied 
publications  for  NSMB  Research  Group  [24-30]. 

Since  the  use  of  PRECAL  is  solely  for  the  determination  of  the 
transfer  functions  (RAOs)  in  (3.22),  it  is  important  to  recognize  that 
any  similar  deterministic  hydrodynamic  analysis  program  may  be  utilized 
for  computing  these  transfer  functions. 

In  the  context  of  PRECAL,  the  pressure  amplitude  per  unit  wave 
amplitude,  PA,  together  with  the  phase  angle,  <j>,  are  given  for  each 
specified  point  for  every  frequency  of  encounter  (u>)  and  wave  heading 
(9).  Then  the  Hpn  (ico,8)  is  calculated  from  the  expression 

lW'",e>  ■  (3.23a) 
or  in  the  expanded  form, 

^pr/io,,e)  =  PA(<O»0)COS(1,(W»9)  +  iPA(u>,6)sin<t>(u),0)  (3.23b) 

Thus,  the  power  spectrum  of  the  random  pressure  loading  may  now  be 
computed  using  Equation  (3.22). 
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3.4  Determination  of  the  Cross  Spectral 
Density  Matrix  of  the  Applied  Loads 


The  power  spectrum  of  the  random  pressure  load  computed  in  the 
preceding  section  is  a  suitable  and  acceptable  input  for  a  finite  element 
program  to  be  used  for  random  response  analysis.  However,  it  cannot  be 
directly  utilized  in  its  raw  form  because  of  the  distributed  nature  of 
this  pressure  load.  A  procedure  is  required  to  convert  this  distributed 
loading  to  the  loads  at  nodal  (global)  degrees  of  freedom  of  the  finite 
element  model.  In  other  words,  a  consistent  finite  element  representa¬ 
tion  of  the  assembled  cross  spectral  density  of  nodal  forces  correspond¬ 
ing  to  the  distributed  random  pressure  field  must  be  determined.  It  is 
only  this  form  of  description  that  can  be  used  in  the  representation  of 
the  discretized  equations  of  motion  (2.21)  and  the  subsequent  solution 
process  discussed  in  the  third  section  of  Chapter  2. 

The  formulation  of  a  consistent  finite  element  representation  of  a 
distributed  random  load  is  quite  involved.  The  most  rigorous  approach  is 
the  development  and  application  of  special  elements  for  the  modelling  and 
discretization  of  these  random  distributed  loads  considered  as  random 
fields.  These  elements,  popularly  referred  to  as  Stochastic  Finite 
Elements  (SFEM)  or  Random  Field  Finite  Elements  have  been  used  quite  a 
bit,  especially  in  applications  Involving  uncertain  (i.e.  stochastic) 
parameters;  see,  for  example,  the  works  of  Contreras  [31],  Liu  et  al. 
[32],  and  Vanmarcke  et  al.  [33].  Stochastic  finite  elements  are  the  most 
suited  for  the  discretization  of  random  fields  because  their  mathematical 
formulation  accounts  for  the  resolution  of  stochastic  inhomogenities  and 
the  extent  of  spatial  correlations.  In  practical  finite  element 
analysis,  the  implication  is  that  two  types  of  elements  or  discretiza¬ 
tions  are  required:  one  for  the  structure  or  the  governing  equations  of 
motions  and  another  for  the  random  fields  present.  Thus,  the  use  of  SFEM 
may  be  computationally  expensive  and  inconvenient.  As  a  result,  there- 
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fore,  less  rigorous  but  acceptably  reliable  methods  of  obtaining  the 
consistent  finite  element  representations  of  distributed  loadings  are 
preferred  for  applications  in  which  only  the  externally  applied  load  is 
stochastic.  These  latter  methods  are  able  to  utilize  conventional  finite 
elements  used  in  deterministic  analysis  for  the  discretization  of  the 
random  pressure  fields  to  a  reasonable  level  of  accuracy. 

In  deterministic  analysis,  if  the  dynamic  surface  pressure  acting  on 
an  element  at  time  t  is  denoted  as  p(s,t)  where  s  is  the  position  on  the 
element,  then  it  is  well  known  that  the  consistent  finite  element 
representation  of  the  nodal  forces  at  time  t,  {Fe(t)j  corresponding  to 
this  distributed  load  is  given  by 


(Fe(t)}  =  /e[N]Tp(s,t)  dAe  (3.24) 

A 

where  Ae  is  the  area  of  the  element  and  [N]  is  the  matrix  of  element 
shape  functions  relating  the  displacements  at  any  point  to  the  element 
nodal  displacements.  Let  CaJi  denote  the  Boolean  matrix  which  relates 
the  nodal  degrees  of  freedom  of  element  i  to  the  nodal  degrees  of  freedom 
of  the  complete  structure.  If  this  pressure  field  acts  over  a  number  of 
elements  NE  of  the  structure,  the  consistent  load  vector  corresponding  to 
the  complete  global  degrees  of  freedom  is  then  given  by  [9] 


NE 


£F(t)}  =  l  [a]T  /  [N]Tp(s,t)  dA. 

i=l  1  Ai  1  "  1 


(3.25) 


where  [N],  is  the  shape  function  for  element  i  and  Ai  is  the  area  of 
the  element. 


Now  let  us  consider  the  case  in  which  p(s,t)  is  random.  Using  the 
definition  of  the  cross  spectral  density  given  in  Equation  (2.14)  of 
Chapter  2,  it  can  be  readily  shown  that  the  cross  spectral  density  of  the 
nodal  forces  in  Equation  (3.25)  is  given  by 
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NE  NE  T 

[SFF(<o)]  =  l  l  [a]!j  /  CN(si)]!spp(si,sj,w)[N(sj)]jdAidA.[a].(3.26) 
1=1  j=l  A-jAj 

where  Spp(s. ,s^ ,w)  is  the  cross  spectral  density  of  the  random  pressure 
field,  and  A. ,Aj  represent  the  surface  area  of  elements  i,j.  Equation 
(3.26)  gives  the  consistent  finite  element  representation  of  the 
assembled  cross  spectral  density  matrix  for  the  random  pressure  load. 

It  can  be  seen  that  the  evaluation  of  Equation  (3.26)  is  computa¬ 
tionally  intensive.  First,  the  double  product  (area  integration)  over 
pairs  of  elements  has  to  be  performed  and  secondly  the  double  summation 
also  must  be  performed.  Considering  that  this  procedure  has  to  be 
performed  for  every  applied  (forcing)  frequency  of  interest  to  the  user, 
the  computation  costs  could  be  phenomenal  1  This  is  why  the  first 
attempts  at  calculating  [SFF(u>)]  for  distributed  surface  pressures  used 
ad-hoc  approximations  to  represent  the  continuous  field  by  concentrated 
transverse  shear  loads  acting  at  the  node  points. 

The  works  of  Olson  [34]  and  Olson  and  Lindberg  [35]  represent  the 
pioneering  efforts  directed  toward  the  calculation  of  the  consistent 
finite  element  representation  of  distributed  random  loads.  These  works 
are  based  on  the  standard  modal  approach  using  the  mode  shapes  obtained 
from  a  finite  element  representation  of  the  continuum  under  study.  A 
polynomial  representation  over  each  finite  element  of  the  excitation 
cross  spectral  is  then  introduced.  This  permits  a  closed-form  evaluation 
of  the  spatial  Integrations  involved  in  the  determination  of  the  force 
cross  spectral  density  matrix.  The  method  was  applied  for  a  four-degree- 
of-freedorn  beam  element  in  [34]  and  a  twenty-degree-of-freedom  triangular 
plate  element  in  [35].  In  [34],  for  example,  it  was  assumed  that 
spp(xi>xj*w)  acting  on  a  uniform  beam  could  be  represented  as: 

Spp(xi,xj,w)  =  ei  +  e  2 x  i  +  e3Xj  +  e4xixJ. 


(3.27) 
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where  the  parameters  e^e*  are  evaluated  in  terms  of  the  values  of  Spp 
for  the  four  node  points  of  the  two  elements  i  and  j,  taking  them  two  at 
a  time.  This  bilinear  approximation  gives 


NE  NE 


Cspp(u)]  -  X.  X  [a]|  [e^fo]  {f0}^  +  ejffJ.ff,}! 

i=j  j=l  1  J  1  J 

+  e3{f0}itf1}]  +  e^f^.ffjT]  [alj 


where 


CM  =  /a[N(x)]T  dx 


CM  =  /^(x)]7  x  dx 


(3.28) 

(3.29a) 

(3.29b) 


are  the  consistent  load  vectors  for  uniform  and  linearly  varying  loads. 
Note  that  e1-e4  are  functions  of  frequency  and  (f0],  {fx}  are  independent 
of  frequency. 


Dey  [36]  employed  a  lumped-load  method  to  convert  the  pressure 
cross-spectral  density  of  the  nodal  forces.  In  this  approach,  it  was 
assumed  that  the  random  nodal  loads  can  be  obtained  by  multiplying  the 
random  forces  by  the  area  allocated  at  the  nodal  points.  In  other  words, 
it  is  assumed  that  the  same  random  force  acts  on  the  entire  area  with 
full  correlation.  This  formulation  gives  the  consistent  finite  element 
cross  spectral  density  matrix  as 


[Spp(u))]  -  £A]  [Spp (s, <o_)]  [A]7 


(3.30) 


in  which  [A]  is  a  diagonal  matrix  of  the  areas  associated  with  each 
global  degree-of-f reedom. 
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In  a  more  recent  study,  Yang  and  Kapania  [37]  presented  a  consistent 
finite  element  formulation  for  a  48  degree-of-freedom  quadrilateral  shell 
element  with  bi-cubic  Hermitian  polynomial  interpolation  functions  as  the 
displacment  shape  functions.  These  shape  functions  are  used  to  form  the 
matrix  of  cross  spectral  densities  of  the  generalized  nodal  forces  for 
distributed  loads  which,  unlike  the  work  of  [36],  are  not  assumed  to  be 
fully  correlated.  The  use  of  Gaussian  quadrature  in  this  work  allows  the 
spectral  density  function  to  be  used  directly  in  its  original  form  rather 
than  in  an  approximate  polynomial  form  as  was  done  in  [34,35].  Another 
interesting  feature  of  the  work  in  [37]  is  that  the  element  shape  func¬ 
tions  are  also  used  to  interpolate  the  spectral  density  values  of  an 
arbitrary  pair  of  points  within  two  separate  shell  elements  from  the 
element  nodal  values.  The  formulation  was  used  to  predict  the  stationary 
random  response  of  cooling  tower  type  shells  in  which  the  random  distri¬ 
buted  loads  were  assumed  as  stationary  in  time  but  permitted  to  be  non- 
homogeneous  in  space.  It  can  be  seen  that  this  formulation  is  actually 
an  implementation  of  Equation  (3.26)  for  the  particular  case  of  the  shell 
element  used  in  their  study. 

3.4.1  Procedure  Implemented  in  this  Work 
for  Random  Load  Description _ 

Rather  than  attempt  to  implement  one  of  the  above  procedures,  a 
different  approach  has  been  devised  during  the  course  of  this  work.  The 
pressure  transfer  functions  are  converted  to  the  nodal  force  transfer 
functions  referred  to  global  degrees-of-freedom  of  the  finite  element 
model.  The  procedure  consists  of  calculating  the  resultant  force  acting 
on  an  element  by  multiplying  the  pressure  by  the  area  of  the  element. 
Then  the  portion  of  this  resultant  force  that  is  allocated  to  each  node 
is  determined  on  the  basis  of  the  consistent  finite  element 
representation  corresponding  to  a  unit  pressure  load  acting  uniformly 
over  the  surface  of  the  element.  The  direction  cosines  of  the  normal  at 
the  element  centroids  are  used  to  compute  the  components  of  the  nodal 
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forces  at  each  node  with  respect  to  the  global  coordinate  axes.  While 
this  procedure  is  short  of  an  exact  consistent  finite  element 
representation  (since  it  will  not  induce  forces  in  the  rotational 
degrees-of-f reedom) ,  it  is  believed  to  be  capable  of  giving  very  good 
approximations,  while  being  more  computationally  efficient.  Indeed,  for 
planar  elements,  the  procedure  yields  exact  results  for  the  three 
elements  implemented,  and  in  any  case  should  always  give  better  results 
than  lumped-load  method  of  Dey  [36]  described  above. 

Once  the  transformation  from  pressure  transfer  functions  to  nodal 
force  transfer  functions  has  been  accomplished,  the  relation: 

SFF(U,)  =  'oX  n(iu’0)HF  n(iw’0)snn<w’0>d0  (i’  j  =  1 . NS) 

1  J'1  11  (3.31) 

is  used  to  directly  compute  the  required  assembled  cross-spectral  density 
matrix  corresponding  to  the  global  degrees-of-f reedom.  It  may  be  noted 
that  the  right  hand  side  of  Equation  (3.31)  is  identical  to  that  of 
Equation  (3.22)  except  that  the  nodal  force  transfer  function  HFn(iu>,0) 
has  replaced  the  pressure  transfer  function  Hpn(ito,0). 

It  is  fortunate  that  this  approach  could  be  devised  and  adopted  in 
this  work  as  the  applications  are  specifically  aimed  at  ship  structures 
in  random  seas.  We  believe  it  is  still  desirable,  however,  to  provide  a 
more  generalized  capability  for  modelling  distributed  random  loads  in 
VAST. 
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4*°  COMPUTER  PROGRAMS  FOR  RANDOM  RESPONSE  ANALYSIS 
4.1  Program  RANVIB 

A  random  vibration  analysis  module  called  RANVIB  has  been  developed 
to  provide  a  random  response  capability  in  the  finite  element  program 
VAST.  The  program  has  been  developed  on  the  basis  of  the  modal  frequency 
response  method  that  was  extensively  discussed  in  Chapter  2.  In  this 
chapter,  a  description  of  the  flow  and  manner  of  operation  of  the  program 
and  its  subprograms  is  presented.  The  main  features  of  the  capabilities 
of  the  program  are  highlighted  and  the  few  limitations  are  noted. 

To  perform  a  finite  element  analysis  using  VAST,  appropriate  Master 
Control  Codes  must  be  selected  and  the  random  vibration  analysis  is  no 
exception.  The  Master  Control  Code  that  triggers  a  random  response  anal¬ 
ysis  is  the  newly  introduced  (eleventh)  code  IPOSTP  which  should  be  set 
equal  to  1,  together  with  a  header  "RANRES"  in  the  PREFX.USE  file.  If 
random  stress  and/or  strain  responses  are  of  interest,  then  additionally 
the  header  ''NODSTR"  should  also  be  present  in  the  PREFX.USE  file.  This 
is  because  the  statistical  properties  of  random  stresses  and  strains  are 
computed  at  the  nodes  using  the  averaged  modal  stresses  or  strains  at  the 
nodal  points  computed  by  P0STV2.  It  should  be  noted  that  all  other  ten 
Master  Control  codes  must  be  selected  to  perform  operations  that  are 
required  before  the  random  response  analysis  is  performed.  For  a  virgin 
VAST  run  for  the  purpose  of  performing  random  response  analysis,  for 
example,  the  control  codes  IELEMS,  IASSEM,  ISTIFM,  IMASSM,  IDECOM, 
IEIGEN,  ILOADS,  IDISPS  should  be  set  to  1.  If  either  random  stresses  or 
strains  are  also  of  interest,  then  the  Master  Control  Code  ISTRES  should 
also  be  set  to  1.  Finally,  IPOSTP  must  be  set  to  1  and  appropriate  head¬ 
ers  provided  in  the  PREFX.USE  file  as  described  above.  Figure  4.1  illus¬ 
trates  the  various  VAST  modules  required  for  a  virgin  random  vibration 
analysis.  Data  preparation  for  the  modules  shown  in  this  diagram  follow. 
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of  course,  exactly  the  procedures  described  in  the  VAST  User's  Manual 
[38].  The  minor  differences  are  discussed  below. 

The  difference  in  data  preparation  starts  with  the  module  DISP5  (see 
Figure  4.1).  For  this  module,  LTYPE  is  set  to  -1  instead  of  1  (for 
steady  state  displacements  due  to  sinusoidal  loading)  or  2  (for  steady 
state  displacements  due  to  a  Fourier  series  representation  of  a  periodic 
loading).  Furthermore,  for  the  damping,  the  modal  damping  (IDAMP  =  1) 
follows  the  same  procedure  as  for  a  frequency  response  analysis.  How¬ 
ever,  the  random  vibration  analysis  also  permits  the  specification  of 
proportional  damping  (IDAMP  =  2)  in  which  only  two  parameters  and 
ctK  (instead  of  NM  modal  damping  ratios)  are  required  to  describe  the 
damping  in  the  manner  described  in  Chapter  2.  For  this  type  of  damping, 
the  user  must  still  specify  NM  values  of  modal  damping  ratios  but  must 
ensure  that  the  first  two  damping  ratios  are  and  respectively, 
so  that 

aM  =  DAMP(l)  (4.1a) 

aK  =  DAMP (2)  (4.1b) 

All  other  modal  damping  ratios  are  ignored  for  this  type  of  damping 
because  they  are  not  needed  in  the  computation  of  the  complex  frequency 
response  functions.  Similarily,  for  structural  damping,  (IDAMP  =  3),  the 
user  inputs  the  structural  damping  factor  (g)  as  DAMP  (1)  while  DAMP(2, 
...  NM)  values  are  entered  as  zero.  All  other  data  for  the  DISP5  module 
is  prepared  as  if  LTYPE=1  as  described  in  the  VAST  user  manual. 

The  load  module  for  a  random  response  analysis  is  different  from 
that  of  a  deterministic  analysis  in  that  the  function  of  the  random  load 
module  is  to  generate  the  consistent  cross  spectral  density  matrix  (and 
not  a  consistent  load  vector)  of  the  nodal  forces  from  the  user  specified 
power  spectra  of  point  loads  and/or  distributed  loads.  A  new  load  module 
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called  RANLOD  has  been  developed  during  the  course  of  this  work  to 
perform  this  function.  A  random  load  data  file  called  PREFX.RLD  which 
describes  the  point  loads  and/or  distributed  loads  for  the  various 
forcing  frequencies  to  be  used  in  the  analysis  must  be  prepared  by  the 
user.  This  data  file  is  then  used  by  RANLOD  to  generate  the  assembled 
consistent  cross  spectral  density  matrix  of  the  applied  loads  which  are 
written  on  a  binary  file  PREFX.PSD  for  each  forcing  frequency.  This  file 
must  be  created  by  running  RANLOD  before  a  random  vibration  analysis  can 
be  performed.  A  description  of  the  random  load  generation  programs  is 
given  in  Section  4.2. 

Next,  if  either  the  random  stresses  or  strains  are  of  interest  to 
the  user,  then  the  modal  stresses  are  computed  and  written  on  the  binary 
file  PREFX.T53.  Modifications  have  been  made  to  the  STRESS  module  to 
permit  the  computation  of  these  modal  stresses.  Modifications  are  also 
to  be  made  to  the  postprocessing  program  P0STV2  to  allow  for  the 
processing  of  the  PREFX.T53  containing  modal  stresses  so  that  averaged 
modal  stresses  and/or  strains  at  the  nodes  can  be  made  available  on 
PREFS.P21  and  PREFS.P22  files. 


Finally,  the  random  vibration  analysis  module  (RANVIB)  is  run  to 
compute  the  random  dynamic  responses  of  interest  to  the  user  as  specified 
in  the  PREFX.RIN  file  which  is  the  input  data  file  for  this  program.  A 
description  of  this  data  file  is  given  in  Appendix  A. 

Program  RANVIB  is  the  main  driver  for  the  random  response  analysis 
subprograms.  It  consists  of  four  main  subprograms:  PRERDM,  MODRES, 
RANRES,  and  PRNTRR.  At  the  beginning  of  a  RANVIB  run,  all  temporary 
files  to  be  used  for  storing  intermediate  results  and  all  permanent  files 
to  be  used  for  saving  the  computations  (and  which  will  be  used  for  graph¬ 
ics  post  processing)  are  opened  and  named. 
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Briefly,  the  functions  of  the  above  named  subprograms  are  as  fol¬ 
lows.  Subroutine  PRERDM  prepares  information  necessary  for  the  computa¬ 
tions  to  be  performed;  subroutine  MODRES  computes  the  random  modal 
response;  subroutine  RANRES  directs  and  controls  the  computation  of  the 
statistical  properties  of  the  quantities  of  interest;  and  subroutine 
PRNTRR  prints  out  the  random  response  results.  A  simplified  flow  diagram 
of  program  RANVIB  is  shown  in  Figure  4.2. 

Subroutine  PRERDM  first  reads  the  PREFX.RIN  input  data  file  for  user 
information.  This  user  information  consists  of  the  physical  processes 
and  statistical  properties  of  interest  to  the  user.  As  described  in 
Appendix  A,  the  first  card  on  the  PREFX.RIN  files  contains  the  parameters 
IRESP  and  IDAMP.  If  IRESP=0,  only  the  modal  response  will  be  computed 
and  saved  on  the  PREFX.RMR  file  which  may  be  used  in  a  future  run.  This 
modal  response  consists  of  the  cross  spectral  density  matrix  of  the  modal 
displacements  [Sqq(w)]  for  all  forcing  frequencies  and  the  covariance 
matrix  of  modal  displacements  [Cqq].  Both  of  these  quantities  are  very 
important  for  random  response  computations.  If  IRESP=1,  the  modal 
response  will  be  freshly  computed  (i.e.  in  the  current  run)  and  will  be 
used  for  random  response  computations.  If  IRESP=2,  a  modal  response 
computed  during  a  previous  run  and  saved  on  PREFX.RMR  file  is  used  for 
the  response  computations.  This  is  like  a  restart  and  saves  a  lot  of 
computation  time.  The  parameter  IDAMP  specifies  the  type  of  damping. 
Other  user  information  on  the  PREFX.RIN  file  includes  the  physical  pro¬ 
cesses  and  statistical  properties  of  interest,  whether  cross  statistics 
or  only  auto-statistics  are  of  interest,  and  whether  responses  for  all 
nodes  or  only  a  selected  number  of  nodes  are  to  be  computed.  Subroutine 
PRERDM  extracts  the  relevant  natural  frequencies  and  orthonormal i zed  mode 
shapes  from  the  PREFX.T51  file  to  form  the  modal  matrix  [R]  which  is 
written  columnwise  on  tape,,  The  PREFX.T52  file  is  read  for  information 
concerning  additional  forcing  frequencies  and  damping. 
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The  random  physical  processes  are  divided  into  two  groups.  The 
first  group  consists  of  displacements,  velocities,  and  accelerations  and 
are  referred  to  as  primary  responses.  The  second  group  which  consists  of 
stresses  and  strains  are  called  secondary  responses.  The  computation  of 
the  statistical  properties  of  secondary  response  processes  has  not  yet 
been  implemented.  If  primary  random  responses  are  needed  only  at  a 
selected  number  of  nodes,  a  reduced  modal  transformation  matrix  [P] 
corresponding  to  these  nodes  is  formed  and  written  columnwise  on  tape 
NTPR1.  If  stresses  are  of  interest,  the  transformation  matrix  [T]  for 
the  stresses  at  the  selected  nodes  or  all  the  nodes  of  the  finite  element 
model  is  determined  from  the  PREFS.P21  file  and  written  columnwise  on 
tape  NTPR2.  Similarly,  the  transformation  matrix  [W]  for  the  strains  is 
determined  from  the  PREFS.P22  file  and  written  columnwise  on  tape  NTPR3. 
It  should  be  noted  that  while  different  statistical  properties  of 
processes  within  a  group  may  be  of  interest,  the  locations  at  which 
properties  of  processes  within  a  group  are  to  be  computed  must  be  the 
same  for  the  processes.  This  is  to  facilitate  the  use  of  the  same  [P] 
for  primary  responses.  For  secondary  responses,  this  allows  the  use  of 
the  same  information  to  be  communicated  to  P0STV2  for  the  generation  of 
PREFS.P21  and  PREFS.P22  files  and  the  possibility  of  generating  both 
files  in  a  single  program  run,  to  be  used  in  the  determination  of  [T]  and 
CVQ«  Furthermore,  this  permits  the  use  of  the  same  number  of  computa¬ 
tional  degrees-of-freedom  (NSPRIM)  for  all  primary  response  processes  and 
the  same  number  of  computational  degrees-  of-freedom  (NSECND)  for  both 
stresses  and  strains.  Another  advantage  of  performing  computations  at 
the  same  locations  for  all  primary  responses  is  that  the  computation  of 
the  spectral  densities  of  one  primary  response  process  from  another  is 
straightforward. 

The  last  operation  performed  by  subroutine  PRERDM  is  the  checking  of 
the  user  requests  to  determine  additional  statistical  properties  that 
will  be  needed  internally  by  the  program  but  have  not  been  explicitly 
requested  by  the  user  in  the  input  file.  For  example,  if  the  user  does 


4.6 


not  request  the  power  spectrum  of  velocities  but  wants  the  correlation 
functions  of  velocities  to  be  computed,  subroutine  PRERDM  will  inform 
RANVIB  that  the  power  spectrum  of  velocities  will  be  computed  by  setting 
IVV(l)  (see  Appendix  A)  equal  to  2  instead  of  the  zero  value  assigned  by 
the  user.  In  this  case,  however,  the  results  of  the  velocity  power 
spectra  will  neither  be  saved  in  the  response  file  nor  printed  in  the 
VAST  output  file  (PREFX.LPT).  The  flow  diagram  of  subroutine  PRERDM  is 
shown  in  Figure  4.3. 

The  next  subprogram  that  is  run  after  PRERDM  is  subroutine  MODRES 
which  performs  the  computation  of  the  modal  responses  according  to  the 
value  of  IRESP  as  described  above.  The  first  operation  within  MODRES  is 
the  arrangement  of  the  forcing  frequencies  in  ascending  order  of  magni¬ 
tude.  This  is  because  random  response  computations  involves  numerical 
integrations  with  respect  to  frequency  and  so  it  is  important  that  the 
frequencies  are  in  order  so  that  the  numerical  algorithm  can  be  applied 
in  order  at  the  different  steps.  It  should  be  noted  that  the  arrangement 
of  the  cross  spectral  density  of  the  loading  on  the  PREFX.PSD  file  is 
also  in  ascending  order  of  magnitude  of  each  forcing  frequency  (including 
the  natural  frequencies  of  interest  to  the  user).  The  natural  frequen¬ 
cies  are  usually  included  In  the  response  computations  so  as  to  be  able 
to  predict  the  peak  responses  accurately  at  these  resonant  frequencies. 
This  is  why  the  natural  frequencies  must  be  known  to  the  user  before 
RANLOD  is  run,  as  is  shown  in  Figure  4.1. 

If  IRESP  is  not  equal  to  2,  MODRES  loops  over  all  forcing  frequen¬ 
cies  to  perform  the  following  computations.  First,  the  cross  spectral 
density  of  the  applied  loads  [Spp(u>)]  for  the  current  frequency  (to)  is 
read  from  the  PREFX.PSD  file.  Then  the  cross  spectral  density  of  the 
modal  force  [Sff(u))]  is  computed  using  Equation  (2.26).  Next,  the 
complex  frequency  response  function  for  the  current  frequency  is  computed 
for  all  the  NM  modes.  This  computation  is  based  on  Equation  (2.29a)  for 
IDAMP=1 ,  Equation  (2.29c)  for  IDAMP=2,  and  Equation  (2.29e)  for  IDAMP=3. 
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Next,  the  cross  spectral  density  matrix  of  the  modal  displacement 
[Sqq(u))]  for  the  current  frequency  is  computed  using  Equation  (2.28), 
and  written  on  tape  NTMODS.  After  the  end  of  the  forcing  frequencies 
loop,  the  covariance  matrix  of  the  modal  displacements  [Cqq]  is 
computed  using  Equation  (2.49)  and  also  written  on  NTMODS  (i.e.  the 
PREFX.RMR  file).  This  function  is  provided  by  subroutine  MODCOV  which  is 
also  capable  of  computing  other  modal  covariances  such  as  modal 
velocities,  modal  accelerations,  modal  acceleration-rates  or  any  other 
modal  rates  through  a  parameter  JFACT.  This  subroutine  has  been  designed 
with  the  objective  that  it  will  be  useful  for  the  computation  of 
higher-order  spectral  moments  (of  any  order)  that  may  be  used  in  spectral 
fatigue  analysis  in  future.  A  simplified  diagram  of  the  flow  of 
subroutine  MODRES  is  illustrated  in  Figure  4.4. 

As  can  be  seen  in  Equation  (2.26)  or  its  indicial  representation  in 
Figure  4.4,  the  computation  of  the  cross  spectral  density  matrix  of  the 
modal  force  from  that  of  the  applied  force  and  the  modal  matrix  is  very 
intensive  since  the  evaluation  of  a  triple  matrix  product  has  to  be 
performed  at  every  forcing  frequency  of  interest.  This  is  even  more  so 
since  neither  the  (NSxNM)  modal  transformation  matrix  [R]  nor  the  NAF 
(NSxNS)  cross  spectral  density  matrices  [Spp(w)]  can  be  stored  in-core 
in  the  computer  for  large  values  of  NM,  NAF,  or  NS  thereby  necessitating 
a  lot  of  input-output  operations.  With  the  exception  of  modal  matrices, 
all  cross  spectral  density  matrices  are  stored  on  tape  as  vectors  (either 
rowwise  or  columnwise)  and  only  two  complex  vectors  are  stored  in-core  at 
any  given  time  during  the  analysis.  Even  for  frequency-dependent  modal 
matrices,  only  one  matrix  at  a  particular  frequency  is  stored  in-core  at 
any  given  time.  This  is  why  the  results  of  the  modal  response  quantities 
(even  though  usually  not  of  interest  to  the  user)  are  saved  on  PREFX.RMR 
so  that  it  can  be  used  in  future  runs  when  additional  statistical  proper¬ 
ties  or  properties  of  other  random  processes  are  required.  Indeed,  in  a 
typical  run,  the  computation  of  the  modal  responses  may  consume  more  CPU 
time  than  the  computation  of  the  statistical  properties  of  the  processes. 


4.8 


An  approximation  may  be  employed  to  reduce  the  computations  involved 
in  the  evaluation  of  [Sff(o))].  This  approximation  assumes  that 
[Sff(u>)]  is  a  diagonal  matrix  whose  diagonal  terms  are  computed  as 

Sfj(w)  =  Wj[SFF(co)]{4>}r  (4.2) 

where  SF£(w)  is  the  diagonal  term  corresponding  to  the  rth  mode,  and  {$} 
is  the  eigenvector  for  mode  r.  However,  this  approximation  is  known  to 
be  valid  only  for  lightly  damped  structures  whose  natural  frequencies  are 
also  widely  separated  [9].  This  approximation  is  not  used  in  RANVIB. 

The  next  subprogram  run  after  MODRES  is  subroutine  RANRES.  This 
subprogram  directs  the  computation  of  the  statistical  properties  of  the 
random  responses  in  the  manner  illustrated  in  Figure  4.5.  First,  the 
modal  displacement  covariance  matrix  [Cqq]  is  read  from  the  PREFX.RMR 
binary  file  on  tape  NTMODS.  Then,  the  spectral  density  of  the  primary 
responses  are  computed  as  required  by  the  user.  This  is  performed  first 
for  displacements,  then  for  velocities  and  lastly  for  accelerations.  The 
other  statistical  properties  of  the  primary  responses  are  then  computed. 
Again,  the  order  of  computation  is:  displacements,  velocities,  accelera¬ 
tions.  The  order  of  computation  of  statistical  properties  is:  covari¬ 
ances,  correlations,  apparent  frequencies.  For  the  secondary  responses, 
stresses  are  treated  first,  then  strains  follow  as  needed.  The  order  of 
computations  of  the  statistics  of  the  responses  in  both  cases  is  spectral 
densities,  covariances,  correlations,  apparent  frequencies.  Finally  the 
statistical  properties  of  nodal  forces  are  computed  following  the  same 
order  as  for  the  secondary  responses.  It  will  be  noted  that  in  all  cases 
the  spectral  densities  are  computed  first  when  needed.  This  is  because 
they  may  be  required  in  the  calculation  of  other  statistical  properties. 
Similarly,  the  covariances  are  computed  before  the  apparent  frequencies 
because  the  former  represent  the  denominator  in  the  definition  of  the 
latter.  Equation  (2.60). 
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Figures  4.6a  and  4.6b  give  the  flow  diagrams  of  the  procedures  for 
computing  the  spectral  densities  of  primary  responses  and  secondary 
responses  respectively.  Both  procedures  form  subroutine  POWERS  which 
proceeds  in  the  form  of  Figure  4.6a  or  4.6b  depending  on  the  value  of  the 
parameter  ICODE.  For  primary  responses,  IC0DE=2,  for  stresses  IC0DE=3 

arnnVn  T'"*'  haS  been  made  for  "°dal  forces 

(IC0DE=1)  which  will  be  implemented  in  the  future.  It  can  be  seen  from 
igures  4.6a  and  4.6b  that  the  major  operation  involved  in  the  computa¬ 
tion  of  the  cross  spectral  densities  is  the  evaluation  of  the  triple 
matrix  products  for  each  of  the  forcing  frequencies,  namely. 


CSxx((o)]  =  CP]CSqq(o))][P]T 


(4.3) 


for  displacements  for  example.  Thus,  computation  time  can  be  saved  if 
the  size  of  [P]  (i.e.  NSPRIM)  is  very  small  and  this  is  accomplished  by 
reques  ing  only  a  selected  number  of  nodes  instead  of  asking  for  respons¬ 
es  at  all  nodes  of  the  model.  The  same  also  applies  to  other  responses 

eluding  stresses  and  strains  when  the  computation  of  secondary 
responses  becomes  operational. 


The  flow  diagrams  illustrating  the  procedures  for  the  computation  of 
covariances  are  shown  in  Figures  4.7a-4.7e.  Different  flow  diagrams  have 
been  provided  for  displacements,  velocities,  accelerations,  stresses,  and 
strains  since  the  procedures  depend  (slightly)  on  the  response  quantity 
o  interest  as  decided  by  the  parameter  ICODE.  A  value  of  IC0DE.1  is  for 

iconr  ,f0rC? '  IC°DE'2  U  f°r  d,Sp,acements'  1C0DE.3  Is  for  velocities, 
IC0DE.4  is  for  accelerations,  ICODE-5  is  for  stresses,  and  ICODE-6  is  for 

s  rains.  All  procedures  are  performed  by  subroutine  COVARC.  An  import¬ 
ant  issue  in  the  computation  of  covariances  is  the  decision  about  whether 
the  covariances  will  be  computed  from  the  corresponding  (or  related) 
spectral  densities  or  from  an  appropriate  modal  covariance  matrix.  A 
was  discussed  in  Chapter  2,  it  is  preferable  to  compute  the  covariance 
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matrix  of  a  response  process  from  the  modal  covariance  matrix  if  the 
spectral  densities  of  the  process  have  not  been  previously  computed 
because  of  the  savings  in  computation  time.  The  triple  matrix  product 
involved  in  this  is  performed  by  subroutine  COVMUL.  However,  if  the 
spectral  densities  are  available,  it  is  better  to  use  them  in  the  compu¬ 
tation  of  the  covariances  because  the  computations  required  are  consider¬ 
ably  reduced  and  they  should  also  give  more  accurate  values  since  the 
spectral  densities  are  directly  integrated.  This  philosophy  is  the  basis 
for  the  flow  diagrams  presented  in  Figures  4.7a-4.7e.  It  will  be  noted 
in  Figures  4.7a-c  that  in  the  cases  concerning  primary  responses,  even  if 
the  spectral  densities  of  the  particular  process  of  interest  is  not 
available  but  those  of  another  primary  response  process  are  available, 
then  the  required  spectral  densities  are  first  readily  computed  from 
those  available  and  then  used  for  determining  the  required  covariances. 
Note  that  in  Figures  4.7a-c,  C*J,  cjjj,  denote  the  elements  of  the 
covariance  matrices  of  modal  displacements,  modal  velocities,  and  modal 
accelerations  respectively.  The  diagonal  terms  of  the  covariance  matrix 
give  the  mean  square  responses  and  the  square  roots  of  those  mean  squares 
gives  the  root  mean  square  responses  commonly  referred  to  as  rms  respons¬ 
es.  In  RANVIB,  the  mean  square  responses  (and  not  rms  values)  are  saved 
on  results  files  and  printed  out  in  the  VAST  output  file. 

The  computation  of  autocorrelation  functions  is  performed  by  subrou¬ 
tine  CORELA  described  in  Figure  4.8.  This  is  essentially  a  straight 

forward  implementation  of  Equations  (2.41),  which  in  numerical  form  has 
the  form: 


Rxx(t) 


i_  NAf_1  fSXX(aW  ~  SXxK)  r 

2irT  i=i  1  («1+1  _  lcos(w1+1x)-cos(w1T)] 


+  SXX{wi+l)s1n(wi+iT)  "  Sxx(a,i)sin(o,iT)}  (4-4) 
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The  sequence  of  NTAU  values  of  t  of  interest  to  the  user  are  indicated 
either  explicitly  or  otherwise  in  the  input  data  file  as  described  in 
Appendix  A.  This  numerical  approximation  assumes  a  linear  variation  of 
SXX(«J)  between  the  discrete  frequencies  o>.  and  w.+r  This  is  the  same 
approximation  assumed  when  the  mean  square  response  is  computed  from  the 
power  spectrum  using: 

1  NAF-1 

E[X2]  =  4tt  tWl)  +  W^Vl  "  ®i>  (4-5) 
as  an  approximation  to  Equation  (2.42). 

Thus,  it  is  important  that  a  sufficient  number  of  frequencies  are  select¬ 
ed  in  the  frequency  range  of  interest  in  order  that  accurate  values  of 
correlation  functions  and  covariances  be  computed  with  the  above  trape¬ 
zoidal  schemes. 

Subroutine  APFREQ  computes  the  apparent  frequencies  as  shown  in 
Figure  4.9.  As  discussed  in  Chapter  2,  the  denominator  in  the  definition 
of  the  apparent  frequency  is  simply  the  mean  square  response  while  the 
numerator  is  the  mean  square  response  of  the  corresponding  rate  process. 
Covariances  required  are  computed  in  the  manner  described  earlier.  Since 
only  autostatistics  are  required  for  apparent  frequencies,  a  special 
subroutine  COVINT  is  available  for  numerically  integrating  autospectra 
(if  previously  computed)  using  Equation  (4.5)  for  determining  mean  square 
responses.  This  subroutine  is  especially  designed  for  use  by  subroutine 
APFREQ  but  is  also  used  by  subroutine  COVARC  if  cross  statistics  are  not 
of  interest  and  autospectra  have  been  previously  computed. 


All  random  response  results  requested  explicitly  by  the  user  are 
saved  in  the  appropriate  files.  In  addition  to  the  PREFX.RMR  file  for 
saving  the  modal  responses,  six  additional  permanent  binary  files  are 
used  to  save  responses.  These  are  PREFX.RFC  for  forces,  PREFX.RDP  for 


displacements,  PREFX.RVL  for  velocities,  PREFX.RAC  for  accelerations, 
PREFX.RSS  for  stresses,  and  PREFX.RSN  for  strains.  Separate  files  are 
used  for  saving  these  results  because  the  volume  of  output  data  written 
on  any  of  the  files  may  be  large,  especially  if  cross  statistics  are 
requested.  These  files  are  later  used  for  plotting. 

The  final  subprogram  run  in  RANVIB  is  subroutine  PRNTRR  which 
directs  the  printout  of  the  random  response  results  in  the  VAST  output 
file  (PREFX.LPT).  This  subroutine  consists  of  subroutines  PRNTPS  for 
printing  the  autospectra,  PRNTMS  for  printing  mean  square  responses, 
PRNTCR  for  printing  autocorrelation  functions,  and  PRNTAF  for  printing 
the  apparent  frequencies.  The  cross  statistical  properties  of  responses 
are  not  printed  out  but  are  available  in  the  permanent  files  to  be  used 
for  post  processing  as  needed.  A  description  of  how  to  request  the 
printouts  of  the  random  response  results  is  included  in  Appendix  A. 
Figure  4.10a  illustrates  the  order  in  which  the  processes  are  printed  out 
while  Figure  4.10b  shows  the  order  in  which  the  response  statistics  of  a 
given  process  are  printed  out. 

4.2  Computer  Programs  for  Random  Load  Generation 

A  set  of  computer  programs  for  generating  the  consistent  finite 
element  cross  spectral  density  matrix  of  the  random  load  for  a  ship  whose 
hydrodynamic  pressure  transfer  functions  are  computed  by  the  HPCFEM  suite 
of  programs  were  developed  to  meet  the  requirements  of  Section  2.0  of 
this  contract.  Figure  4.11  shows  the  skeletal  flow  diagram  of  the 
procedure  for  generating  the  PREFX.PSD  file  which  is  essentially  the 
random  load  information  required  by  program  RANVIB  for  random  response 
analysis. 


The  three  programs  developed  during  the  course  of  this  work  are 
referred  to  as  PREHPC,  POSHPC,  and  RANLOD. 
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Program  PREHPC  is  a  pre-interface  program  to  the  HPCFEM  suite  of 
programs.  The  main  function  of  this  program  is  to  utilize  information 
from  the  VAST  input  data  files  and  its  own  input  data  file  called 
PREFX.HPC  for  the  generation  of  the  input  data  files  required  to  run  the 
HPCFEM  programs.  These  input  data  files  include  HYDMES_SHIP_INPUT.DAT, 
HYDMES_SHIP_GEOMETRY_INPUT.DATA,  PRECAL_INPUT.DAT,  and  FINMES_INPUT.DAT. 
The  procedure  for  creating  the  input  data  file  (PREFX.HPC)  for  running 
PREHPC  is  described  in  Appendix  B.  Subprogram  HGLWD  in  module  HVAST  may 
be  used  for  generating  a  PREFX.LWD  which  contains  the  sectional  mass 
distribution  along  the  ship  longitudinal  direction.  The  procedure  for 
running  program  HVAST  can  be  found  in  [39].  Similarly,  subprogram  DIGHLL 
in  module  HLLFLO  may  be  used  for  generating  the  offset  table  for  the  ship 
lines  of  form  if  the  offsets  of  the  two-dimensional  sections  of  the  ship 
is  not  otherwise  known.  Program  DIGHLL  uses  a  digitizing  tablet  to 
create  an  offset  table  from  an  offset  diagram  and  the  offset  table  is 
stored  in  a  file  PREFX.DGH  for  use  by  program  PREHPC.  Of  course,  there 
are  user  options  to  manually  provide  the  sectional  masses  and  the  offsets 
in  the  input  data  file  PREFX.HPC  if  they  are  known.  Further  details  on 
the  procedure  and  requirements  for  running  DIGHLL  are  given  in  [40]. 

Program  PREHPC  also  generates  an  output  file  PREFX.PRE  which 
contains  information  about  the  ocean  environment  in  which  the  ship  is 
travelling.  Such  information  includes  the  type  of  ocean  (as  identified 
by  a  suitable  ocean  wave  spectrum),  the  number  of  encounter  frequencies, 
the  number  of  wetted  elements,  the  element  type,  the  element  connectivi¬ 
ties,  and  the  direction  cosines  of  the  normals  of  the  elements. 

Next,  the  HPCFEM  suite  of  programs  are  run  in  the  order  HYDMES, 
PRECAL,  and  FINMES.  The  most  relevant  output  file  is  the  file  called 
FINMES_PRESSURE.OUT  which  contains  the  pressure  transfer  functions  at 
points  corresponding  to  the  centroids  of  the  elements. 
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The  next  program  to  be  run  is  POSHPC  which  is  the  post-interface 
program  between  the  HPCFEM  and  VAST  programs.  Program  POSHPC  extracts 
the  pressure  transfer  functions  contained  in  file  FINMES_PRESSURE.OUT  and 
utilizes  this  information  together  with  the  information  on  PREFX.PRE  file 
for  automatically  generating  a  PREFX.LOD  file.  This  PREFX.LOD  file  is 
like  a  load  file  containing  point  loads  for  use  in  a  transient  dynamic 
analysis.  The  number  of  time  steps,  NTIME  is  set  equal  to  (NAF*NWH*2)  in 
which  NAF  is  the  number  of  forcing  (or  encouter)  frequencies  and  NWH  is 
the  number  of  wave  headings.  The  factor  2  accounts  for  the  fact  that  the 
pressure  transfer  function  at  any  point  is  a  complex  number  with  non-zero 
real  and  imaginary  parts  so  that  the  phase  information  is  preserved.  The 
procedure  for  generating  the  point  loads  used  to  create  the  PREFX.LOD 
file  is  an  approximate  procedure  for  computing  the  consistent  nodal 
forces  corresponding  to  a  given  distributed  loading  by  using  a  nodal 
load-lumping  technique  derived  from  the  knowledge  of  the  exact  consistent 
finite  element  load  vector.  In  fact,  the  approximate  procedure  is  exact 
for  planar  elements  and  was  implemented  for  three  elements  in  the  VAST 
element  library.  These  elements  include  the  thick-thin  shell  element 
(IEC=1) ,  the  triangular  plate  element  (IEC=4),  and  the  quadrilateral 
shell  element  (IEC=5).  In  essence,  therefore,  the  pressure  transfer 
functions  are  converted  into  nodal  force  transfer  functions  corresponding 
to  the  global  degrees  of  freedom  of  the  finite  element  model.  The 
PREFX.LOD  file  when  used  to  run  program  L0AD1  module  in  VAST  generates  a 
binary  load  file  (PREFX.T47)  which  contain  the  nodal  force  transfer 
functions  to  be  later  used  by  program  RANLOD. 

The  subprogram  RANLOD  is  responsible  for  the  computation  of  the 
cross  spectral  density  matrix  corresponding  to  the  given  random  loading. 
RANLOD  utilizes  information  on  the  PREFX.PRE  and  PREFX.T47  files  for 
computing  the  complex  nodal  force  transfer  functions  for  all  the  wave 
headings  and  forcing  frequencies  concerned.  The  cross  spectral  density 
matrix  is  then  computed  at  each  forcing  frequency  (o>)  by  using  the 
expression: 


4.15 


sff(°))  -  ■^o1THFin^ia,,0^HF.n^iO5,e^snn^lo,G^d0  (1»>1»NS) 

J 

where  Spjl(w)  denotes  the  cross  spectral  density  of  the  force  for  degrees- 
of-freedom  i  and  j,  Hp^Ow.G)  is  the  complex  nodal  force  transfer 
function  for  i  at  wave  heading  (0),  H?.n(ico,0),  and  Snn(»,e)  is  the 
value  of  the  ocean  wave  spectrum  evaluated  at  the  current  frequency  and 

wave  direction.  The  skeletal  flow  diagram  for  RANLOD  is  shown  in  Figure 
4.12. 


The  program  currently  permits  the  use  of  nine  types  of  ocean  wave 
spectra  which  characterize  the  type  of  ocean  in  which  the  ship  is  travel¬ 
ling  or  stationed.  These  ocean  wave  spectra  (identified  by  the  parameter 
ISPEC)  include  the  following: 

I SPEC  =  1  :  TEN-PARAMETER  JUSZKO  SPECTRUM 

=  2  :  SIX-PARAMETER  OCHI-HUBBLE  SPECTRUM 
=  3  :  P I ERSON-MOSKOWI TZ  SPECTRUM  (W.R.T.  WIND  SPEED) 

=  4  :  PIERSON-MOSKOWITZ  SPECTRUM  (W.R.T.  PEAK  FREQUENCY) 

=  5  :  INTER.  SHIP  STRUCTURES  CONGRESS  (ISSC)  SPECTRUM 
=  6  :  JONSWAP  SPECTRUM  (W.R.T.  GRAVITATIONAL  ACCLN.) 

=  7  :  JONSWAP  SPECTRUM  (W.R.T.  SIGNIFICANT  WAVE  HEIGHT) 

=  8  :  INTER.  TOWING  TANK  CONFERENCE  (ITTC)  SPECTRUM 
-  9  :  BRETSCHNEIDER  SPECTRUM 

Descriptions  of  these  ocean  wave  spectra  were  given  in  Section  3.2 
of  Chapter  3.  A  summary  of  their  mathematical  expressions  and  the 
parameters  required  for  their  description  is  given  in  Appendix  B,  where 
the  Input  data  file  PREFX.HPC  is  described. 


FIGURE  4.1.  A  Flow  Diagram  Showing  Various  VAST  Modules 

Required  for  a  Virgin  Random  Vibration  Analysis 


FIGURE  4.2.  Simplified  Flow  Diagram  of  Program  RANVIB 


FIGURE  4.3.  Flow  Diagram  of  Subroutine  PRERDM  used  for  Preparing 
Information  for  Random  Vibration  Analysis 


FIGURE  4.4.  Flow  Diagram  of  Subroutine  MODRES  for 
Computing  the  Cross  Spectral  Densities 
and  Covariance  Matrix  of  Modal  Displacements 


Note:  Primary  responses  are  computed  in  the  order:  DISPLACEMENTS, 
VELOCITIES,  ACCELERATIONS 

Statistics  are  computed  in  the  order:  PSD,  COV,  CORR. ,  AP.  FREQ. 


FIGURE  4.5.  Flow  Diagram  of  Subroutine  RANRES  which  Directs  the 
Computation  of  the  Response  Statistics  for  all 
Processes 


FIGURE  4.6a.  Flow  Diagram  for  the  Computation  of  the 

Spectral  Densities  of  Primary  Responses  in 
Subroutine  POWERS 


I  I 


FIGURE  4.6b.  Flow  Diagram  for  the  Computation  of  the 
Spectral  Densities  of  Secondary  Responses 
in  Subroutine  POWERS 


Figure  4.7a.  Flow  Diagram  for  the  Computation  of  the 
Covariance  Matrix  of  Displacements  in 
Subroutine  COVARC 


START 


Figure  4.7b.  Flow  Diagram  for  the  Computation  of  the  Covariance 
Matrix  of  Velocities  in  Subroutine  COVARC 


Figure  4.7c.  Flow  Diagram  for  the  Computation  of  the  Covariance 
Matrix  of  Accelerations  in  Subroutine  COVARC 


Figure  4.7d.  Flow  Diagram  for  the  Computation  of  the  Covariance 
Matrix  of  Stresses  in  Subroutine  COVARC 


Figure  4.7e.  Flow  Diagram  for  the  Computation  of  the  Covariance 
Matrix  of  Strains  in  Subroutine  COVARC 


START 


FIGURE  4.8.  Flow  Diagram  of  Subroutine  CORELA  for  the 
Computation  of  Correlation  Functions 


FIGURE  4.9.  Flow  Diagram  of  Subroutine  APFREQ  for  the 
Computation  of  Apparent  Frequencies 


FIGURE  4.10a.  Flow  Diagram  of  Subroutine  PRMTRR  for  Printing  the 
Response  Results  in  the  VAST  Output  File  (PREFX.LPT) 


START 


FIGURE 


.10b.  Flow  Diagram  Illustrating  the  Order  in  which 
the  Various  Response  Statistics  are  Printed 
out  for  Each  Process 


.  Procedure  for  Generating  the  Consistent  Cross  Spectral 
Density  Matrix  of  the  FEM  of  a  Ship  whose  Hydrodynamic 
Pressure  Transfer  Functions  are  Computed  by  the  HPCFEM 
Suite  of  Programs 


FIGURE  4.11 


START 


RUN  PRECAL/FINMES  TO  GENERATE  NODAL 
PRESSURE  TRANSFER  FUNCTIONS 
(RAOs  FOR  VARIOUS  w's  and  6‘s) 


DO  FOR  APPLIED  FREQUENCIES  (w's) 
1  TO  NAF 


DO  FOR  WAVE  HEADING  (0's) 
1  TO  NWH 


RUN  LOAD1  FOR  REAL  PART  OF  l.e.  Use  nodal  pressures 


NODAL  FORCES  RAOs  (HFn(u,0)) 


PA(w,e)cosit>(w,0) 
to  obtain  H*®(w,B) 
by  running  LOAD1 


RUN  LOAD1  FOR  IMAGINARY  PART  OF  l.e.  use  nodal  pressures 
NODAL  FORCE  RAOs  P.(w,0) s1n<J>(u),0) 


(l.e.  H*JJ  (u,0) 


to  obtain  H^'(u,0) 
by  running  LOAD1 


FORM  COMPLEX  NODAL  FORCE  RAOs  FROM  ABOVE  RESULTS: 

Hf  n(1w,0)  -  H^en(u,0)  +  1Hjmn(u),0) 

t  k  k 

(k  =  1,2 . NS) 


END  OF  WAVE  HEADING  LOOP 


GENERATE  X-PSD  OF  NODAL  FORCES  FOR  THE 
CURRENT  FREQUENCY  USING 


;FF<W>  "  '0  Hpin(1«,e)HFjn(1U,0)Snn(W,0)d0 
(1,j  =  1,2 . NS) 


WRITE  CURRENT  [SF(-(w)]  ON  PREFX.PSD 


END  OF  APPLIED  FREQUENCIES  LOOP 


FIGURE  4.12.  Procedure  for  Generating  the  PREFX.PSD  File 
for  Random  Wave  Loads 


5.1 


5.0  GRAPHICS  SUPPORT  FOR  RANDOM  VIBRATION  ANALYSIS 

This  aspect  of  the  contract  has  been  deferred  to  the  next  fiscal 

year. 


The  graphics  support  capabilities  that  are  planned  are  essentially 
"x-y"  plots  of  input  and  output  quantities.  Such  plots  include  the 
power  spectral  density  of  pressures,  forces,  displacements,  velocities, 
accelerations,  stresses,  strains,  versus  frequency  on  linear, 
logarithmic,  or  semi-logarithmic  scales.  It  is  also  intended  to  plot 
correlation  functions  of  output  quantities  as  functions  of  time  delay. 

Furthermore,  in  order  to  be  able  to  convey  a  quick  overall  picture 
of  the  extent  of  deformation,  it  is  planned  to  provide  the  capability  of 
plotting  stress  and  strain  contours  using  the  appropriate  RMS  values  of 
the  random  responses. 


6.1 


6.0  EXAMPLE  PROBLEMS 

A  simple  example  problem  for  the  verification  of  the  random 
vibration  analysis  program  is  the  five  span  continuous  beam  exposed  to 
jet  noise  excitation  as  shown  in  Figure  6.1.  This  problem  has  been 
analyzed  by  several  workers  including  Olson  [34].  This  example  was 
selected  because  it  provides  a  rigorous  verification  of  the  computer 
program.  First,  there  are  non-zero  off-diagonal  terms  in  the 
cross-spectral  density  matrix  of  the  applied  load  arising  from  the 
distributed  nature  of  the  random  pressure  field  induced  by  the  jet  noise 
propagation.  Secondly,  for  the  five-bay  beam,  the  natural  frequencies 
are  not  widely  separated  as  can  be  seen  in  Table  6.1.  Thus,  the 
inclusion  of  the  cross-modal  term  contributions  to  the  overall  response 
is  important. 

The  five  span  straight  beam  is  simply  supported  at  its  ends  and  at 
the  four  intermediate  supports  as  shown  in  Figure  6.1.  Each  span  has 
unit  length.  The  beam  has  a  bending  stiffness  of  unity,  a  mass  per  unit 
length  of  10-4  units,  and  a  cross-sectional  area  of  1000.0  units.  Each 
span  is  modelled  with  four  general  beam  elements  (IEC=3).  The  damping  is 
taken  as  1%  of  critical  damping  for  each  of  the  15  natural  modes  of 
vibration  used  for  computing  the  random  response. 

The  jet  noise  excitation  popularly  referred  to  as  plane  wave 
propagation  of  clipped  white  noise  is  assumed  to  propagate  along  the  beam 
in  the  direction  shown  in  Figure  6.1.  Its  power  spectral  density  Is: 

S(<o,X)  =  exp(-iwX/C0L)  (6.i) 

where  X  =  X2-Xj  is  the  distance  between  the  points  of  interest  located  at 
Xx  and  X2,  C0  is  the  nondimensional  speed  of  propagation  taken  to  be  6 
(following  Olson),  and  L  is  the  length  of  a  span  of  the  beam  which  is 
taken  to  be  1. 


6.2 


The  power  spectral  densities  of  displacements  obtained  in  the  VAST 
run  agree  very  well  with  the  results  of  Olson  £34],  ABAQUS  [11]  and 
Johnsen  and  Dey  [11.36].  For  nodes  3,  7,  11,  15,  and  19,  the  power 
spectral  density  of  displacements  are  given  in  Table  6.2.  Just  as  in 
ABAQUS  and  the  paper  by  Johnsen  and  Dey,  the  VAST  results  in  Table  6.2 
were  computed  by  a  procedure  in  which  the  equivalent  nodal  forces  corres¬ 
ponding  to  the  distributed  loading  were  approximated  by  lumping  at  the 
nodes.  Although  this  procedure  gives  acceptably  accurate  final  results 
for  the  power  spectral  densities  and  mean  square  responses,  an  examina¬ 
tion  of  some  intermediate  modal  force  response  results  show  that  the  use 
of  a  consistent  finite  element  representation  of  the  distributed  random 
load  gives  more  accurate  results.  This  is  because  the  consistent  finite 
element  representation  ensures  the  inclusion  of  terms  in  the  cross- 
spectral  density  matrix  induced  by  loading  at  rotational  degrees-of- 
freedom  which  are  not  accounted  for  in  a  load  lumping  approximation. 
Table  6.3  gives  a  comparison  of  the  results  of  the  first  modal  forces  for 
the  first  ten  forcing  frequencings  considered  in  this  analysis.  The 
exact  results  were  given  by  Olson  [34]  as: 


S1 1 
ff 


(w)  = 


2it2A2 


[ir2-(u>/C0)2] 


[1  +  cos(5w/C0)] 


(6.2) 


where  A  is  the  exact  normalized  amplitude  of  the  first  mode  which  Is 

0.63246  for  this  beam.  It  can  be  seen  that  the  percentage  errors 

associated  with  the  lumping  load  procedure  are  considerably  more  than 

those  obtained  in  a  VAST  analysis  in  which  a  consistent  finite  elelment 

representation  of  the  cross-spectral  density  of  the  random  loads  is 

used.  It  can  also  be  noted  that  the  errors  increase  as  the  values  of  the 

forcing  frequencies  increase  because  the  exact  mathematical  expression 

given  in  (6.2)  is  very  sensitive  to  the  value  of  ».  This  demonstrates 

the  importance  of  the  capability  to  model  distributed  random  loads  a 

capability  which  is  lacking  in  ABAQUS  and  the  ASKA  II  program  used  by  Dey 

C36J. 


6.3 


The  mean  square  responses  for  the  translational  displacement 
degree-of-freedom  and  the  rotational  degree-of-freedom  (i.e.  slope)  are 
shown  in  Tables  6.4A  and  6.4B,  respectively.  It  can  be  seen  that  the 
VAST  analysis  results  agree  very  well  with  the  results  of  Olson  [34].  It 
is  also  interesting  to  note  that  the  mean  square  responses  obtained  in 
VAST  using  the  two  different  methods  described  in  Chapter  4  gave  exactly 
the  same  results.  In  other  words,  the  mean  square  responses  computed 
directly  from  the  power  spectral  densities  and  those  computed  from  the 
modal  mean  square  response  gave  the  same  results.  This  is  a  further 
check  on  the  program. 

The  tight  budget  did  not  permit  the  possibility  of  testing  more 
example  problems.  However,  a  simple  two-degree-of-freedom  system  given 
on  Page  361  of  the  book  by  Elishakoff  [7]  was  the  first  problem  used  for 
debugging  and  testing  the  program  and  the  correct  answers  were  obtained. 
It  is  hoped  that  more  example  problems  will  be  run  in  the  future. 


TABLE  6.1 


NATURAL  FREQUENCIES  OF  FIVE  BAY  BEAM 


Mode  No 


Natural 

Frequency  (rad/s) 

Exact 

Computed  by  VAST 

9.8696 

9.8671 

10.9498 

10.9478 

13.6927 

13.6917 

17.2469 

17.2505 

20.764 

20.7169 

39.478 

39.5564 

41.731 

41.8316 

46.903 

47.0642 

53.124 

53.3782 

58.936 

59.3007 

88.826 

90.0506 

92.182 

93.5441 

99.798 

101.4734 

108.795 

110.7851 

117.056 

119.1292 

TABLE  6.2 

POWER  SPECTRAL  DENSITIES  OF  DISPLACEMENTS  (in  di stance2/HZ) 


Node 

FREQ 

=  1.57 

HZ 

FREQ 

=  1.74 

HZ 

FREQ 

=  2.18 

HZ 

No. 

Johnsen 

& 

Dey 

ABAQUS 

VAST 

Johnsen 

& 

Dey 

ABAQUS 

VAST 

Johnsen 

& 

Dey 

ABAQUS 

VAST 

3 

0.0989 

0.0963 

0.0968 

0.8447 

0.7911 

0.7917 

0.1127 

0.0845 

0.0848 

n 

0.0848 

0.0879 

0.0879 

0.3251 

0.3419 

0.3412 

0.0170 

0.0229 

0.0231 

ii 

0.0609 

0.0607 

0.0610 

0.0018 

0.0021 

0.0022 

0.1696 

0.1701 

0.1709 

15 

0.0411 

0.0375 

0.0378 

0.3115 

0.3002 

0.2991 

0.0170 

0.0110 

0.0110 

19 


0.0318  0.0342  0.0344 


0.8222  0.8882  0.8859 


0.1127  0.1409  0.1414 


TABLE  6.3 


POWER  SPECTRAL  DENSITY  OF  MODAL  FORCES 
FOR  THE  FIRST  MODE  (SFM(1,1))  AT  VARIOUS  FORCING  FREQUENCIES 


Forcing 

Frequency 

(rad/s) 

Exact 

Consistent 

FE  Load 

Representation 

Load  Lumping 
Approximation 

Value 

% 

Error 

Value 

% 

Error 

0.0 

0. 162114E+00 

0. 161999E+00 

0.07 

0.145764E+00 

10.08 

0.986711E+01 

0.983786E-01 

0.983152E-01 

0.06 

0.911085E-01 

7.39 

0. 109478E+02 

0.833137E-02 

0.833969E-02 

0.10 

0.778139E-02 

6.60 

0. 136917E+02 

0.509409E+00 

0.508663E+00 

0.15 

0.484238E+00 

4.94 

0. 172505E+02 

0.234608E+01 

0.234470E+01 

0.06 

0.230600E+01 

1.71 

0.207169E+02 

0. 184701E+01 

0. 184906E+01 

0.11 

0. 188878E+01 

2.26 

0.395564E+02 

0.715720E-02 

0.759556E-02 

6.12 

0. 103953E-01 

45.24 

0.418316E+02 

0.238297E-03 

0.257319E-03 

7.98 

0.36696E-03 

53.99 

0.470642E+02 

0.310576E-02 

0.354043E-02 

14.00 

0.555000E-02 

78.74 

0.533782E+02 

0.308945E-02 

0.388650E-02 

25.80 

0.678970E-02 

119.77 

TABLE  6.4A 


ROOT  MEAN  SQUARE  DISPLACEMENT  RESPONSE 


Node  No. 

RMS  Displacement 

Olson  [34] 

VAST 

%  Difference 

1 

0.0 

0.0 

0.00 

2 

0.1719 

0.1702 

3 

0.2274 

0.2257 

4 

0.1557 

0.1535 

1.41 

5 

0.0 

0.0 

0.0 

6 

0.1225 

0.1197 

2.29 

7 

0.1534 

0.1517 

1.11 

8 

0.1040 

0.1029 

1.06 

9 

0.0 

0.0 

0.0 

10 

0.0904 

0.0944 

4.42 

11 

0.1176 

0.1246 

5.95 

12 

0.0841 

0.0879 

4.52 

13 

0.0 

0.0 

0.0 

14 

0.0889 

0.0878 

1.24 

15 

0.1360 

0.1341 

1.40 

16 

0.1129 

0.1107 

1.95 

17 

0.0 

0.0 

0.0 

18 

0.1585 

0.1581 

0.25 

19 

0.2391 

0.2399 

0.33 

20 

0.1793 

0.1800 

0.39 

21 

0.0 

0.0 

0.0 

TABLE  6.4B 


ROOT  MEAN  SQUARE  SLOPE  RESPONSE 


Node  No 


rm: 

5  Displacement 

Olson  [34] 

VAST 

%  Difference 

0.7988 

0.7894 

1.18 

0.5101 

0.5062 

0.76 

0.2775 

0.2628 

5.30 

0.5319 

0.5300 

0.36 

0.6308 

0.6116 

3.05 

0.3436 

0.3409 

0.79 

0.2421 

0.2280 

5.82 

0.3662 

0.3613 

1.34 

0.4378 

0.4388 

0.23 

0.2819 

0.2986 

5.92 

0.2253 

0.2241 

0.53 

0.2902 

0.3064 

5.58 

0.3801 

0.3827 

0.68 

0.3308  1 

0.3265 

1.30 

0.2216 

0.2143 

3.29 

0.3005 

0.2972 

1.10 

0.6113 

0.6004 

1.78 

0.5652 

0.5684 

0.57 

0.2198 

0.2159 

1.77 

0.5378 

0.5397 

0.35 

0.8235 

0.8263 

0.34 

Jet  Noise  Propagation 
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FIGURE  6.1.  A  Five  Bay  Beam  Excited  by  Jet  Noise 


7.1 


7-0  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 


The  integrity  checks  and  verification  problems  in  respect  of  the 
random  vibration  analysis  program  has  been  very  encouraging.  Overall,  a 
quite  respectable  random  vibration  analysis  capability  has  been  provided 
for  the  VAST  finite  element  program.  The  capability  that  has  been 
provided  surpasses  that  of  most  commercial  finite  element  programs  in 
that  cross  statistical  properties  are  computed,  a  specialized  distributed 
random  load  modelling  capability  has  been  implemented,  and  a  variety  of 
statistical  properties  of  different  processes  may  be  requested  by  the 
user. 


There  are  a  few  items  that  remain  to  be  implemented  to  complete  the 
first  phase  of  the  contract.  First,  the  procedure  for  computing  the 
statistical  properties  of  stress  and  strain  responses  could  not  be 
implemented  due  to  budget  limitations.  However,  an  elegant  procedure  has 
been  designed  to  achieve  this  end  in  the  near  future.  Secondly,  the 
testing  of  the  random  load  generation  program  using  a  real  ship  structure 
could  not  be  performed  because  there  was  no  provision  made  for  that  in 
the  contract.  Finally,  the  provision  of  a  graphics  support  capability 
has  been  deferred  to  the  next  fiscal  year. 

In  addition  to  the  above,  it  is  recommended  that  immediate 
enhancements  to  the  RANVIB  module  should  include  the  following: 

(i)  Capability  for  analysis  involving  the  use  of  superelements; 

(11)  Capability  for  response  analysis  of  structures  excited  by 
support  motion  so  that  problems  involving  earthquake 
excitations  (for  example)  can  be  dealt  with; 

(iii)  Capability  for  dealing  with  models  that  make  use  of  multipoint 
constraints; 


7.2 


(iv)  Provision  of  a  more  generalized  capability  for  computing  the 
consistent  cross  spectral  density  matrix  of  distributed  random 
loads  as  discussed  in  Chapter  3;  and 

(v)  General  algorithmic  improvements  to  increase  computation 
speeds. 

On  the  long  term,  the  development  of  RANVIB  should  be  vigorously 
pursued  because  probabilistic  analysis  is  gaining  wider  acceptance  among 
designers  and  analysts.  In  the  future  other  capabilities  to  be  provided 
should  include  the  ability  to  handle  combined  deterministic  and  random 
loading,  situations  involving  initial  stresses,  and  response  to 
nonstationary  and  nonhomogeneous  loads. 

It  should  be  noted  that  the  main  purpose  of  a  probabilistic  analysis 
of  a  structure  is  to  be  able  to  assess  its  integrity  and,  above  all,  make 
quantitative  statements  about  the  probability  that  it  will  perform  its 
intended  function  satisfactorily.  The  subject  that  deals  with  this  is 
referred  to  as  structural  reliability  and  so  the  ultimate  objective 
should  be  the  provision  of  a  reliability  analysis  capability  in  VAST.  A 
reliability  analysis  module  would  make  use  of  the  response  statistics 
computed  in  the  random  vibration  analysis  module  to  make  quantitative 
predictions  about  the  probability  of  failure  or  safety  of  the  structure. 
The  scientific  authority,  recognizing  that  the  issue  of  fatigue  damage 
and  reliability  are  very  important  to  ship  structures,  included  a 
provision  in  the  second  phase  of  the  contract  for  conducting  a  state-of- 
the-art  literature  review  and  making  appropriate  recommendations.  This 
aspect  of  the  contract  was  also  performed  during  this  year  and  a  separate 
report  [41]  on  the  literature  review  and  recommendations  has  been 
prepared  as  called  for  in  the  contract. 

It  is  hoped  that  the  scientific  authority  maintains  an  interest  in 
this  project  and  continues  to  give  it  the  dedicated  commitment  and 
support  it  truly  deserves. 
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A.  1 


APPENDIX  A 

DESCRIPTION  OF  INPUT  DATA  FILE  PREPARATION 
_ FOR  RUNNING  PROGRAM  RANVIB 


The  following  is  a  description  of  the  method  for  creating  the 
PREFX.RIN  file  to  be  used  as  the  input  data  file  for  running  a  random 
response  analysis  in  VAST. 

The  user  is  advised  to  exercise  caution  in  requesting  the 
computations  of  the  various  statistical  properties  of  the  different 
possible  random  processes  because  a  random  response  analysis  is  more 
expensive  than  a  corresponding  deterministic  dynamic  analysis.  As  such, 
it  is  recommended  that  only  responses  of  interest  should  be  selected  by 
the  user  through  the  appropriate  codes.  Also,  where  possible,  selective 
computations  in  which  only  a  selected  number  of  nodal  locations  are 
requested  are  recommended. 


FORMAT  OF  INPUT  FILE:  PREFX.RIN 


A. 2 


NOTE:  STATS  =  Statistical  Properties 

Primary  Responses  =  displacements,  velocities,  and  accelerations 

Secondary  Responses  =  stresses  and  strains 

(Computation  of  secondary  responses  not  yet  implemented.) 


Card  1A  (215) 

I  OPT 

=  i 

Modal  frequency  response  method  for  analysis 

=  2 

Direct  frequency  response  method  for  analysis 

I  DAMP 

=  1 

Modal  damping  ratios  used  for  describing  damping 
(€y,  y=l ,2, . . . ,NM) 

(NOTE:  Not  applicable  for  I0PT=2) 

=  2 

Proportional  or  Rayleigh  Damping  (C  =  OmM-m^K) 

=  3 

Structural  Damping  (CX  =  igKX) 

Card  IB 

(15)  (Omit 

if  I0PT*1) 

IRESF1 

=  0 

Only  modal  response  to  be  computed  during  this  run, 
i.e.  for  the  sole  purpose  of  creating  PREFX.RMR 

=  1 

Random  response  to  be  computed  using  the  modal 
response  computed  during  the  current  program  run 

=  2 

Random  response  to  be  computed  using  the  modal 
response  computed  during  a  previous  program  run 
(PREFX.RMR  must  be  available  for  this  selection. 
This  is  like  a  restart  and  will  usually  reduce 
computation  time  significantly.) 

Card  1C  (15)  (Omit  if  IOPT/2) 

IRESP2  =  0  Only  complex  frequency  response  function  to  be 

computed  during  this  run,  i.e.  for  the  sole  purpose 
of  creating  PREFX.CFR 

=  1  Random  response  to  be  computed  using  the  complex 

frequency  response  function  computed  during  the 
current  program  run 

=  2  Random  response  to  be  computed  using  the  complex 

frequency  response  function  computed  during  a 
previous  program  run  (PREFX.CFR  must  be  available 
for  this  selection.  This  is  like  a  restart  and 
will  usually  reduce  computation  time  significant¬ 
ly.) 


Card  ID  (2E10.3)  (Omit  if  IDAMP  *  2) 


ALFAM  =  a**  (constant  in  Rayleigh  damping  model) 
ALFAK  =  a|<  (constant  in  Rayleigh  damping  model) 

Card  IE  (E10.3)  (Omit  If  IDAMP  *  3) 

STRDMP  =  g  (structural  damping  factor) 


Card  IF  (415) 
I FORCE  =  0 
=  1 
=  -1 

IPRIMY  =  0 
=  1 
=  -1 
ISECND  =  0 
=  1 
=  -1 
NTAU  =  0 

=  NTAU 


Force  STATS  not  required 

Force  STATS  wanted  at  all  nodes 

Force  STATS  wanted  only  at  selected  nodes 

STATS  of  primary  responses  not  required 

STATS  of  primary  responses  wanted  at  all  nodes 

STATS  of  primary  responses  wanted  at  selected  nodes 

STATS  of  secondary  responses  not  required 

STATS  of  secondary  responses  wanted  at  all  nodes 

STATS  of  primary  responses  wanted  at  selected  nodes 

No  correlation  functions  are  required 

Correlation  functions  will  be  required  in  the 
analysis  (maximum  of  100) 

where 

NTAU  =  number  of  sequence  of  non-zero  t  values  for 
which  the  correlation  functions  are  to  be 
computed 

(Note:  If  ISECND*0,  the  user  must  have  requested 
the  computation  of  modal  stresses  by  set¬ 
ting  the  master  control  code  ISTRES  to  1  so 
that  the  PREFX.T53  file  containing  modal 
stresses  would  have  been  generated  and  used 
in  P0STV2  for  the  generation  of  averaged 
nodal  stresses/strains.  Of  course,  the 
user  must  also  have  appropriately  requested 
a  P0STV2  batch  run.) 


C.§rd  2  Information  for  the  calculation  and  printout  of  the  statistical 
properties  (i.e.  STATS)  of  random  nodal  forces 

(Omit  if  IFORCE=0) 


Card  2A  (15) 

IFX  =  1  Only  AUTO-STATS  of  forces  are  of  interest 

=  2  Both  AUTO-  and  CROSS-STATS  of  forces  are  of  inter¬ 

est 


Card  2B  (415) 
IFST(l)  =  0 
=  1 

=  -1 

IFST (2)  =  0 
=  1 
=  -1 
IFST(3)  =  0 
=  1 

=  -1 

IFST (4)  =  0 
a  1 

=  -1 


Power  spectral  densities  of  forces  not  of  interest 

Power  spectral  densities  of  forces  to  be  saved  on 
random  force  file  and  printed  in  output  file 

Power  spectral  densities  of  forces  saved  but  not 
printed 

Covariance  of  forces  not  of  Interest 

Covariance  of  forces  to  be  computed  and  printed 

Covariance  of  forces  to  be  computed  but  not  printed 

Autocorrelation  functions  of  forces  not  of  interest 

Autocorrelation  functions  of  forces  to  be  computed 
and  printed 

Autocorrelation  functions  of  forces  to  be  computed 
but  not  printed 

Apparent  frequencies  of  force  not  of  interest 

Apparent  frequencies  of  force  to  be  computed  and 
printed 

Apparent  frequencies  of  force  to  be  computed  but 
not  printed 


Card  2C  (omit  if  IFORCE*-!) 


NFN 


Number  of  nodes  for  which  nodal  force  STATS 
required 


are 


(JFN(I) ,  1=1,  NFN) 


Actual  node  numbers  (arranged  In  ascending  order) 
for  which  STATS  of  nodal  forces  are  wanted. 

(Note:  Error  message  is  issued  if  IF0RCE*0  but 
IFST(1)=IFST (2)=IFST (3)=IFST (4)=0) 


-ar(*  -  Information  for  the  calculation  and  printout  of  the  STATS  of 
primary  responses 

(omit  if  IPRIMY=0) 


Card  3A  (415) 
IPX  =  0 

=  1 

I  PD  =  0 

=  1 

IPV  =  0 

=  1 

IPA  =  0 

=  1 


Only  AUTOSTATS  of  primary  responses  are  of  interest 

Both  AUTO-  and  CROSS-STATS  of  primary  responses 
wanted 

Displacement  STATS  not  wanted 
Displacement  STATS  wanted 
Velocity  STATS  not  wanted 
Velocity  STATS  wanted 
Acceleration  STATS  not  wanted 
Acceleration  STATS  wanted 

(Note:  Error  message  is  issued  if  IPRIMY*0  but 
IPD=IPV=IPA=0) 


Card  3B(i)  (415)  (omit  if  IPD=0) 

IDD(1)  =  0  P°wer  spectral  densities  of  displacements  not  of 

interest 

=  1  Power  spectral  densities  of  displacements  to  be 

computed  and  printed 


A. 6 


=  -1 

IDD(2)  =  0 
=  1 

=  -1 

IDD(3)  =  0 
=  1 

=  -1 

IDD(4)  =  0 
=  1 

=  -1 


Power  spectral  densities  of  displacements  to  be 
computed  but  not  printed 

Covariance  of  displacements  not  of  interest 

Covariance  of  displacements  to  be  computed  and 
printed 

Covariance  of  displacements  to  be  computed  but  not 
printed 

Correlations  of  displacements  not  of  interest 

Correlations  of  displacements  to  be  computed  and 
printed 

Correlations  of  displacements  to  be  computed  but 
not  printed 

Apparent  frequencies  of  displacements  not  of  inter¬ 
est 

Apparent  frequencies  of  displacements  to  be  comput¬ 
ed  and  printed 

Apparent  frequencies  of  displacements  to  be  comput¬ 
ed  but  not  printed 

(Note:  Error  message  is  issued  if  IPD=1  but 
IDD(1)=IDD(2)=IDD(3)=IDD(4)=0) 


Card  3B(ii 1  (415)  (omit  if  IPV=0) 


IVV(l) 

=  0 

Power  spectral  densities  of  velocities 
interest 

not  of 

=  1 

Power  spectral  densities  of  velocities  to 
puted  and  printed 

be  com¬ 

=  -1 

Power  spectral  densities  of  velocities  to 
puted  but  not  printed 

be  com- 

IVV(2) 

=  0 

Covariance  of  velocities  not  of  interest 

=  1 

Covariance  of  velocities  to  be  computed  and  printed 

-  -1 

Covariance  of  velocities  to  be  computed 
printed 

but  not 

A. 7 


I VV  (3)  = 


IVV (4)  = 


Card  3B  f  i  i 
IAA(l)  = 


IM(2)  = 

=  1 

IAA(3)  =  0 
=  1 

IAA(4)  =  0 
=  1 


Correlations  of  velocities  not  of  interest 

Correlations  of  velocities  to  be  computed  and 
printed 

Correlations  of  velocities  to  be  computed  but  not 
printed 

Apparent  frequencies  of  velocities  not  of  interest 

Apparent  frequencies  of  velocities  to  be  computed 
and  printed 

Apparent  frequencies  of  velocities  to  be  computed 
but  not  printed 

(Note:  Error  message  is  issued  if  IPV=1  but 
IVV(1)=IVV(2)=IVV(3)=IVV(4)=0) 


1  (415)  (omit  if  IPA=0) 

Power  spectral  densities  of  accelerations  not  of 
interest 


1 


1 


Power  spectral  densities 
computed  and  printed 

Power  spectral  densities 
computed  but  not  printed 


of  accelerations  to  be 
of  accelerations  to  be 


Covariance  of  accelerations  not  of  interest 

Covariance  of  accelerations  to  be  computed  and 
printed 

Covariance  of  accelerations  to  be  computed  but  not 
printed 


Correlations  of  accelerations  not  of  interest 

Correlations  of  accelerations  to  be  computed  and 
printed 

Correlations  of  accelerations  to  be  computed  but 
not  printed 


Apparent  frequencies  of  accelerations  not  of  inter- 


Apparent  frequencies  of  accelerations  to  be  comput¬ 
ed  and  printed 


=  -1  Apparent  frequencies  of  accelerations  to  be  comput¬ 

ed  but  not  printed 

(Note::  Error  message  is  issued  if  IPA=1  but 
IAA(1)=IAA(2)=IAA(3)=IAA(4)=0) 

Card  3C  (omit  if  IPRIMY*1) 

Card  3C(i)  (15) 

NDNP  =  Number  of  selected  Displacement  Nodes  for  which 

Primary  responses  are  to  be  computed 

Card  3C f i i 1  (1615) 

JPN  =  Actual  displacement  node  numbers  (arranged  in 

ascending  order)  for  which  primary  responses  are  to 
be  computed 

(JPN(I), 1=1, NDNP) 

Card  4  (omit  if  (ISECND=0) 

Card  4A  (215) 

ISX  =  0  Only  AUTOSTATS  of  secondary  responses  are  of  inter¬ 

est 

=  1  Both  AUTO-  and  CROSS-STATS  of  secondary  responses 

are  of  interest 

ISC  =  1  Stress  STATS  only  are  wanted 

=  2  Strain  STATS  only  are  wanted 

=  3  Both  stress  and  strain  STATS  are  wanted 

Card  4BM1  (415)  (omit  if  ISC=2) 

ISS(l)  =  0  Power  spectral  densities  of  stresses  not  of  inter¬ 

est 

=  1  Power  spectral  densities  of  stresses  to  be  computed 

and  printed 


-1 


Power  spectral  densities  of  stresses  to  be  computed 
but  not  printed 


ISS(2)  = 


ISS (3)  = 


ISS (4)  = 


Card  4B(ii 
ISN(l)  = 


ISN(2)  = 


ISN(3)  = 


ISN(4)  = 
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Covariance  of  stresses  not  of  interest 

Covariance  of  stresses  to  be  computed  and  printed 

Covariance  of  stresses  to  be  computed  but  not 
printed 

Correlations  of  stresses  not  of  interest 

Correlations  of  stresses  to  be  computed  and  printed 

Correlations  of  stresses  to  be  computed  but  not 
printed 

Apparent  frequencies  of  stresses  not  of  interest 

Apparent  frequencies  of  stresses  to  be  computed  and 
printed 

Apparent  frequencies  of  stresses  to  be  computed  but 
not  printed 


1  (415)  (omit  if  ISC=1) 

0  Power  spectral  densities  of  strains  not  of  interest 

1  Power  spectral  densities  of  strains  to  be  computed 

and  printed 

Power  spectral  densities  of  strains  to  be  computed 
but  not  printed 

0  Covariance  of  strains  not  of  interest 

1  Covariance  of  strains  to  be  computed  and  printed 

-1  Covariance  of  strains  to  be  computed  but  not  print¬ 

ed 

3  Correlations  of  strains  not  of  interest 

L  Correlations  of  strains  to  be  computed  and  printed 

_1  Correlations  of  strains  to  be  computed  but  not 

printed 

)  Apparent  frequencies  of  strains  not  of  interest 

Apparent  frequencies  of  strains  to  be  computed  and 
printed 
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Card  4C 
NGNS 

JSEN 


Card  5 

Card  5A 
I COREL 

Card  5B 
RTAU1 

DELTAU 

Card  5C 
RTAU 


=  -1  Apparent  frequencies  of  strains  to  be  computed  but 

not  printed 


(omit  if  ISECND*-!) 

Number  of  selected  geometric  nodes  for  which 
secondary  responses  are  to  be  computed 

Actual  geometric  node  numbers  (arranged  in  ascend¬ 
ing  order)  for  which  secondary  responses  are  to  be 
computed 

(JSEM(I) ,I=1,NGNS) 

(Note:  The  elements  corresponding  to  these  nodes 
must  have  been  communicated  to  the  stress 
module  as  described  in  the  VAST  manual.) 


(omit  if  NTAU=0) 


(15) 

»  1  1st  value  of  x  (non-zero),  RTAU1,  and  constant 

increment  in  values  of  t,  DELTAU,  are  provided 

=  2  Non-zero  NTAU  values  of  t  are  provided 


2E10.3  (omit  if  IC0REL=2) 

1st  (non-zero)  value  of  time  delay  for  the  computa¬ 
tion  of  correlation  functions 

Constant  increment  in  values  of  x  to  be  used  for 
generating  other  values  of  time  delay 


8E10.3  (omit  If  IC0REL=1) 

The  NTAU  values  of  time  delay  for  which  correlation 
functions  are  to  be  computed 
(RTAU(I) ,  1=1,  NTAU) 
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APPENDIX  B 


Input  Data  for  Program  PREHPC 


Card  1A  (A80)  Job  title  for  HPCFEM  suite  of  programs  (maximum  of  80 
characters). 


Card  IB  (A80)  Ship  name  (maximum  of  80  characters) 


Card  2  (15)  (Defines  the  type  of  the  ocean  wave  spectrum,  Snn) 


ISPEC 


=  1  :  TEN-PARAMETER  JUSZKO  SPECTRUM 

=  2  :  SIX-PARAMETER  OCHI-HUBBLE  SPECTRUM 

=  3  :  PIERSON-MOSKOWITZ  SPECTRUM  (W.R.T.  WIND  SPEED) 

=  4  :  PIERSON-MOSKOWITZ  SPECTRUM  (W.R.T.  PEAK  FREQUENCY) 
=  5  :  INTER.  SHIP  STRUCTURES  CONGRESS  (ISSC)  SPECTRUM 
=  6  :  JONSWAP  SPECTRUM  (W.R.T.  GRAVITATIONAL  ACCLN. ) 

=  7  J  JONSWAP  SPECTRUM  (W.R.T.  SIGNIFICANT  WAVE  HEIGHT) 

=  8  :  INTER.  TOWING  TANK  CONFERENCE  (ITTC)  SPECTRUM 
=  9  :  BRETSCHNEIDER  SPECTRUM 


Card  2.1A  (5F13.6)  (Omit  if  ISPEC  *  1) 


PARAM(l)  =  RL(1)  =  Ai 
P ARAM (2)  =  RL(2)  =  A2 
P ARAM (3)  =  WM(1)  =  ow 
P  ARAM  (4)  =  WM(2)  =  uw 
P ARAM (5)  =  DEL(l)  =  S» 
PARAM(6)  =  DEL(2)  =  S2 


Card  2. IB  (4F13.6)  (Omit  If  ISPEC  /  1) 

PARAM(7)  =  P(l)  =  Pl 
P ARAM (8)  =  P(2)  =  p2 
P ARAM (9)  =  THM(l)  =  6mi 
PARAM(IO)  =  THM(2)  =  0m2 

where  PARAM(I) ,  1=1,  2,  ...,  10,  are  the  ten  parameters  that  describe  the 
Juszko  spectrum  whose  model  description  is  given  by  equation  (3  21) 
namely:  '  ’ 


Snri(w»0) 


14*1+1  Whii'J^Si2  e 
r(Ai)w  4Ai+1 


_riAi±i 

1  4 


* 


cos2pi(0-0m.) 

2 


(B.l) 


B.2 


Card  2.2  (6F13.6)  (Omit  if  ISPEC  #  2) 


This  card 
spectrum: 


snn(w)  = 


defines  the  six  parameters  that  define 


%i4)AiSi2 
4  1 _ 

r(Ai)w  4Ai+1 


the  Ochi-Hubble 


(B.2) 


PARAM(l)  =  RL(1)  = 

P ARAM (2)  =  RL(2)  =  A2 
PARAM(3)  =  WM  ( 1 )  =  (0,,,! 
P  ARAM  (4)  =  WM(2)  = 

P ARAM (5)  =  DEL(l)  =  S, 
PARAM(6)  =  DEL(2)  =  S2 


Card  2.3  (3F13.6)  (Omit  if  ISPEC  *  3) 

This  card  defines  the  three  parameters  of  the  Pierson-Moskowitz  spectrum 
described  in  terms  of  the  wind  speed: 

Snn(w)  =  ag2orB  exp  [-0.74(^)-4]  (B  3) 

PARAM(l)  =  ALPHA  =  a  (Fixed  at  0.0081) 

PARAM(2)  =  GAC  *  g  =  gravitational  acceleration 
PARAM(3)  =  VW  =  Uw  =  wind  speed 


Card  2.4  (3F13.6)  (Omit  if  ISPEC  *  4) 

This  card  defines  the  three  parameters  of  the  Pierson-Moskowitz  spectrum 
described  in  terms  of  the  frequency  of  spectral  peak  (<o0): 

Snn(w)  =  ag2or5  exp  [-1.25(^-)-4]  (B>4) 

PARAM(l)  =  ALPHA  =  a  (Fixed  at  0.0081) 

PARAM(2)  ==  GAC  =  g  =  gravitcitional  acceleration 
PARAM(3)  =  WO  =  w0  =  frequency  of  spectral  peak 


Card  2.5  (2F13.6)  (Omit  if  ISPEC  *  5) 

This  card  defines  the  two  parameters  of  the  International  Ship  Structures 
Congress  (ISSC): 


t 


Snn<®>  =  0.1107  H*  gj  exp [-0. 4427(g)-]  , 
where  u>  =  1.296  w0 

PARAM(l)  =  HS  =  Hs  =  significant  wave  height 
PARAM(2)  =  WO  =  W„  =  peak  frequency 


Card  2.6  (4F13.6)  (Omit  if  ISPEC  *  6) 

This  card  defines  the  four  parameters  of  the  JONSWAP  spectrum 
with  respect  to  gravitational  acceleration: 

Snn(w)  =  ag2urB  exp[-1.25  y  exP[“|pw^  3 

UJ0 


PARAM(l)  =  ALPHA  =  a  (Usually  taken  to  be  0.0081) 
PARAM(2)  =  GAC  =  g  =  gravitational  acceleration 
PARAM(3)  =  WO  =  w0  =  peak  frequency 
PARAM(4)  =  GAMA  =  y  =  peakedness  parameter 


Card  2.7  (3F13.6)  (Omit  if  ISPEC  #  7) 

This  card  defines  the  three  parameters  fo  the  JONSWAP  spectrum 
with  respect  to  significant  wave  height: 


Srm(w) 


exp[-1.25(g-)— ] 


Y  exP 


» 


where 


a*  _  _ 0.0624 _ 

0.230  +  0.0336y  -  0.185  (1.9+y)-1 

PARAM(l)  =  HS  =  Hs  =  significant  wave  height 
PARAM(2)  =  WO  =  (o0  =  peak  frequency 
PARAM(3)  =  GAMA  =  y  =  peakedness  parameter 
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(B.5) 


described 


(B.6) 


described 


(B.7a) 

(B.7b) 


Card  2.8  (3F13.6)  (Omit  if  ISPEC  *  8) 

This  card  defines  the  three  parameters  of  the  International  Towing  Tank 
Conference  (ITTC): 


* 


B.4 


S(<o)  =  ag2w-B  expt-^-55— ] 

H  2 


where  a  == 


0.0081 

k“ 


and  k 


To 

3.54w. 


in  which  a  =  i/m^  =  |  H$  . 


(B.8a) 

(B.8b) 

(B.8c) 

(B.8d) 


In  the  above,  o  is  the  standard  deviation  (r.m.s.  value)  of  the  water 
surface  elevation,  mn  is  the  nth  moment  of  the  spectrum  defined  as 


mn  =  [  ®%n<®)<l®  .  (B.8e) 

and  the  average  zero  crossing  frequency  or  apparent  frequency  (defined  in 
Chapter  2),  w2  is: 

®z  =  {  i2-  }%  (B.8f ) 

o 

PARAM(l)  =  ALPHA  =  a  (calculated  as  given  by  (B . 8b) 

PARAM(2)  =  GAC  =  g  =  gravitational  acceleration 
PARAM(3)  =  HS  =  Hs  =  significant  wave  height 

Card  2.9  (2F13.6)  (Omit  If  ISPEC  *  9) 

This  card  defines  the  two  parameters  of  the  Bretscheider  spectrum: 


Snn(u>)  =  0.1687  H l  exp[-0.675(jjj)-4 ] 


(B.9a) 


where  ws  = 


22 

Ts 


Ts  =  0.946To. 


(B.9b) 


PARAM(l)  =  HS  =  Hs  =  significant  wave  height 
PARAM(2)  =  TO  =  T0  =  peak  period. 
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Card  3  (2F15.5;  2E10.3) 

CDEN  =  Density  of  the  seawater  (eg.  1.025  tonnes/m3) 

CACC  =  Acceleration  due  to  gravity  (eg.  9.81  m/s2) 

(User  must  ensure  that  the  units  of  CDEN  and  CACC  are 
consistent  with  units  used  for  defining  the  VAST  finite  element 
model.) 


Card  4  (5F15.5;  5E10.3) 

RLOA  =  Length  overall 

RLBP  =  Length  between  perpendiculars 

BM  =  Maximum  beam 

TD  =  Mean  draught 

TR  =  Trim  (positive  by  the  stern) 

Card  5  (215) 

ISYM  =  0  :  ship  has  no  symmetry 

=  1  :  ship  has  symmetry  about  X=0  only 

=  2  :  ship  has  symmetry  about  Y=0  only 

=  3  :  ship  has  symmetry  about  both  X=0  and  Y=0 

ISPED  =  1  :  ship  speed  is  in  m/s 

=  2  :  ship  speed  is  in  knots 


Card  6  (215) 

IMAS  =  1  :  sectional  mass  data  is  provided  in  this  input  data  file 

=  2  :  sectional  mass  data  computed  from  the  finite  element  mesh 
of  the  ship  is  available  on  file  PREFX.LWD  (generated  by 
subprogram  HGLWD)  of  module  HVAST 

NSEC  =  No.  of  stations  along  the  global  X-axis  at  which  sectional 
masses  are  defined 


£-arc*  Z  (4F15.5;  4E10.3)  (To  be  provided  NSEC  times) (Omit  if  IMAS  =  2) 
SMASS  =  Sectional  mass 

XM  =  X  coordinate  of  the  C.O.G.  of  the  sectional  mass 

YM  =  Y  coordinate  of  the  C.O.G.  of  the  sectional  mass 

ZM  =  Z  coordinate  of  the  C.O.G.  of  the  sectional  mass 

(C.O.G.  =  centre  of  gravity) 
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Card  8  (415)  (Omit  if  ISYM  =  0) 

NAFT1  =  Number  of  nodes  used  to  define  the  aftend  of  the  ship 
NAFT2  =  Number  of  facets  that  make  up  the  aftend  of  the  ship 
NF0R1  =  Number  of  nodes  used  to  define  the  forend  of  the  ship 
NF0R2  =  Number  of  facets  that  make  up  the  forend  of  the  ship 


Card  9  (15,  3E10. 3) (Provide  NAFTI  Cards) (Omit  if  ISYM  =  0) 

NODEN  =  Number  of  identifier  of  the  aftend  node  point 

XCO  =  X  coordinates  of  the  node  point 

YCO  ==  Y  coordinates  of  the  node  point 

ZCO  =  Z  coordinates  of  the  node  point 


Card  10  (515) (Provide  NAFT2  Cards) (Omit  if  ISYM  =  0) 

NF  ==  Facet  number  for  this  aft  end  facet 

(NN1,  NN2,  NN3,  NN4)  =  Connectivities  of  this  facet 

Facets  Cards  (Applicable  to  both  aftend  and  forend  cards) 

NF  is  the  facet  identifier.  1-36  are  the  facet  identifiers  for 
the  facets  on  the  water-line,  starting  from  the  stern  centre-line 
and  rotating  in  a  clockwise  fashion,  looking  from  above.  The 
facets  around  the  girth  (starting  from  the  port-side)  and 
directly  under  the  water-line  facet  with  identifier  N  are 
numbered  N01,  N02,  N03,  ..  etc.  Nl,  N2,  N3,  N4  are  the  node 
numbers  that  make  up  that  facet.  The  numbering  system  is  in 
sequential  clock-wise  order  looking  from  outside  the  ship  and 
starting  from  the  top  right-hand  corner  of  the  facet. 

Card  11  (15,  3E10.3)  (Omit  if  ISYM  =  0) 

NODEN  =  Number  identifier  of  the  forend  node  point 

XCO  ==  X  coordinates  of  the  node  point 

YCO  =  Y  coordinates  of  the  node  point 

ZCO  ==  Z  coordinates  of  the  node  point 


Card  12  (515)  (Omit  if  ISYM  =  0) 

NF  =  Facet  number  for  this  forend  facet  (NN1,  NN2,  NN3,  NN4)  = 
connectivities  of  this  facet. 


Card  13A  (15)  (Omit  if  ISYM  =  0) 

NLONG  =  No.  of  longitudinal  positions  at  which  automatic  panels  are  to 
be  generated  by  HPCFEM  program 


A 
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Card  13B  (3E10.3,  15)  (Omit  if  ISYM  =  0) 

XX  =  Longitudinal  position  (expressed  as  fraction  of  the  length) 

(Note:  0.0  corresponds  to  aft  perpendicular  and  10.0 

corresponds  to  fore  perpendicular) 

NP  =  No.  of  panels  required  to  model  the  half  girth 


Card  14A  (215) 

NSEC2  =  No.  of  two-dimensional  sections  for  which  body  plan  offsets  are 
defined 

IDATA  =  1  :  horizontal  and  vertical  body  plan  offsets  are  provided  in 
this  input  data  file 

2  :  the  horizontal  and  vertical  offsets  are  provided  by  the 
output  file  of  the  digitizing  program  DIGHLL 


Card  14B  (Omit  if  IDATA  =  2)  (To  be  provided  NSEC2  times) 

RSEC(I)  =  Station  number  (Station  No.  0.0  corresponds  to  A.P.;  Station 
No.  10.0  corresponds  to  F.P.) 

Y(I)  =  Horizontal  offset,  measured  from  the  centre-line  starting  at 
the  keel,  all  values  positive 

Z(I)  =  Vertical  offset,  measured  from  the  base-line  starting  at  the 
keel,  positive  upwards 


Card  15  (E10.3)(i.e.  free  format) 
DAMP  =  Damping  ratio  factor 


Card  16  (15,  4E10.3) 

NSPEED  =  No.  of  ship  speeds  (max.  of  4) 
SPED(I)  =  Values  of  the  ship  speeds 
(1=1,  NSPEED) 


Card  17A  (15) 

NHEAD  =  No.  of  ship  headings  (max.  of  13) 


Card  17B  (8E10.3)  (Provide  1  or  2  cards  as  required) 
HEAD  =  Values  of  ship  headings  in  degrees 
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Card  18A  (15) 

NFREQ  =  No.  of  wave  frequencies  (max.  of  15) 

Card  18B  (8E10.3)  (Provide  1  or  2  cards  as  required) 

FREQ(I)  =  Wave  frequencies  in  rad/s 
(1=1,  NFREQ) 

Card  19A  (15) 

NLPFOM  =  0  forces  and  moments  on  ship  not  desired 

0  no  of  longitudinal  positions  for  which  forces  and  moments  on 
ship  are  to  be  calculated  and  output 


Card  19B  (8G10.3)  (Provide  as  many  cards  as  required) 

(XL ( I ) ,  I  =  1,  NLPFOM) 

=  Longitudinal  positions,  expressed  as  a  fraction  of  the 
length,  for  calculation  of  forces  and  moments 
(Note:  A.P.  =  0.0,  F.P.  =  1.0) 


Card  20  (15) 

INTOPT  =  Interpolation  option  parameter 

=  1  the  closest  hydrodynamic  pressure  point  to  required  point 
=  2  weighted  average  of  hydrodynamic  pressure  points  closest  to 
the  required  point 

=  3  linear  interpolation  between  the  hydrodynamic  pressure  points 
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